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PROJECT  SUMMARY 


The  quantitative  prediction  of  the  dynamics  of  separating  unsteady  flows,  such  as 
dynamic  stall,  is  of  crucial  importance.  This  six-month  SBIR  Phase  I  study  has 
developed  several  new  pressure-based  methodologies  for  solving  3D  Navier-Stokes 
equations  in  both  stationary  and  moving  (body-comforting)  coordinates.  The 
present  pressure-based  algorithm  is  equally  efficient  for  low  speed  incompressible 
flows  and  high  speed  compressible  flows.  The  discretization  of  convective  terms  by 
the  presently  developed  high-order  TVD  schemes  requires  no  artificial  dissipation 
and  can  properly  resolve  the  concentrated  vortices  in  the  wing-body  with  minimum 
numerical  diffusion.  It  is  demonstrated  that  the  proposed  Newton's  iteration 
technique  not  only  increases  the  convergence  rate  but  also  strongly  couples  the 
iteration  between  pressure  and  velocities.  The  proposed  hyperbolization  of  the 
pressure  correction  equation  is  shown  to  increase  the  solver's  efficiency.  The  above 
proposed  methodologies  were  implemented  in  an  existing  CFD  code,  REFLEQS. 
The  modified  code  was  used  to  simulate  both  static  and  dynamic  stalls  on  two-  and 
three-dimensional  wing-body  configurations.  Three-dimensional  effect  and  flow 
physics  are  discussed.  Further  development  and  validation  are  proposed  for  Phase 
II. 
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1. 


INTRODUCTION 


1.1  Dynamic  Stall  Phenomenon  and  Its  Significance 

Dynamic  stall  is  a  complex  physical  event  induced  by  a  large  amplitude  motion  of 
aerodynamic  bodies.  It  is  a  phenomenon  characterized  by  the  shedding  and  passage 
over  the  upper  surface  of  a  lifting  surface  of  vortex-like  disturbances.  Associated 
with  this  phenomenon  is  the  generation  of  intense  vorticity  near  the  nose  of  the 
body,  which  occurs  as  the  pitching  of  the  lifting  surface  dynamically  surpasses  its 
stall  angle  of  attack.  This  vorticity  increases  the  circulation  of  the  flow  and  thus  the 
lift  force  acting  on  the  body.  As  a  result,  large  unsteady  aerodynamic  forces  are 
generated  from  which  the  lift,  drag  and  moment  coefficients  greatly  exceed  their 
maximum  static  counterparts.  The  unsteady  effects  of  dynamic  stall  are  usually 
dominated  by  turbulent  flow  and  the  production  of  large  scale  vortices.  The 
dynamic  stall  events  will  either  proceed  with  the  generation  of  weaker  vortices  if 
the  body  remains  pitched  above  its  static  stall  angle  of  attack,  or  terminate  if  the  body 
returns  to  an  angle  of  attack  sufficiently  small  for  reattachment  of  the  flow.  Figure 
1-1  illustrates  the  typical  flow  field  during  dynamic  stall.  Excellent  reviews  on  the 
subject  have  been  presented  by  McCroskey^'^  and  Carr.^ 

Dynamic  stall  is  of  importance  in  various  aerodynamic  applications  including 
aircraft  maneuverability,  helicopter  rotors,  and  wind  turbine.  For  example,  when 
the  dynamic  stall  appears  in  the  retreating  blade  of  a  helicopter  rotor,  it  produces  a 
loss  of  lift,  thus  an  increase  in  power  is  required  which  in  turn  increases  the 
pitching  loads  and  vibratory  stresses.  Therefore,  significant  efforts  have  been 
devoted  to  understand  and  eliminate  the  undesirable  effects  associated  with 
dynamic  stall  on  helicopter  rotors.  Recent  efforts  are  exploring  the  possibility  of 
utilizing  the  unsteadiness  of  the  flow  field  to  enhance  aircraft  performance  and  to 
attain  the  sustained  dynamic  maneuvering  in  the  post-stall  flight  regime.  For 
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example,  Robinson  and  Luttges'^*®  have  analyzed  the  repetitive  interaction  of  the 
dynamic  stall  vortices  as  a  means  of  maintaining  flow  attachment  to  increase  lift  at 
high  angles  of  attack.  They  also  correlated  this  phenomenon  as  a  function  of  the 
driving  parameters  involved,  i.e.,  pivot  location,  airfoil  shape,  Reynolds  number 
and  mean  angle  of  attack.  Indeed,  the  impetus  to  exploit  the  energetic  nature  of 
large  vortices  to  potentially  enhance  performance  has  already  been  demonstrated.^ 
However,  it  is  clear  that  before  such  a  realistic  usage  is  possible,  extensive  studies 
must  be  undertaken  to  expand  our  knowledge  of  fundamental  aspects. 

1.2  Literature  Review 

Dynamic  stall  is  much  more  difficult  to  analyze  and  predict  than  static  stall  because 
it  depends  on  many  parameters,  including; 

•  airfoil  pitch  rate,  pitch  amplitude,  and  pitch  axis  location; 

•  mean  and  maximum  angles  of  ramp  or  oscillation; 

•  airfoil  geometry,  including  thickness,  leading  edge  curvature,  and  camber; 

•  wing  or  blade  tip  shape; 

•  free  stream  Reynolds  and  Mach  numbers. 

In  the  past,  dynamic  stall  research  has  proceeded  along  several  avenues:  analytical, 
experimental,  and  computational.  Analytical  methods  were  used  primarily  to 
complement  both  the  experimental  and  computational  methods.  Due  to  the 
complexity  of  the  flow  field,  analytical  methods  were  often  found  not  self-sufficient, 
and  as  such  will  not  be  discussed  herein. 

As  this  study  concerns  computational  investigation  of  dynamic  stall  phenomenon, 
the  related  works  on  computational  approach,  experimental  visualizations,  and 
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three-dimensional  effect  have  been  reviewed,  and  briefly  summarized  in  the 
following  three  sections. 


1.2.1  Computational  Approach 

The  recent  advancement  of  Computational  Fluid  Dynamics  (CFD)  has  provided  a 
prediction  capability  that  was  unattainable  only  a  few  years  ago.  The  earliest 
computational  investigations  of  dynamic  stall  appeared  in  the  1970s  and  early  1980s 
with  the  works  of  Mehta,^  Gulcat,^  and  Sankar®'^®  as  the  most  representative  ones. 
In  1977  Mehta^  solved  the  Navier-Stokes  equations  for  laminar  flows  to  determine 
the  flow  field  around  an  oscillatory  NACA0012  airfoil.  Although  the  flow  field  was 
computed  in  detail,  it  was  limited  to  low  Reynolds  number  flows  and  required  large 
amounts  of  computer  time. 

Relatively  recently,  Sankar  and  his  co-workers®"^®  used  the  unsteady,  compressible, 
Reynolds  averaged  Navier-Stokes  equations  in  the  computation  of  laminar  and 
turbulent  flow  fields  around  oscillating  airfoils.  Their  aerodynamic  load  prediction 
agreed  with  the  experimental  data  during  the  upstroke  prior  to  the  stall.  Yet,  they 
could  not  resolve  the  details  in  the  post-stall  regime. 

In  the  mid  1980s,  Wu  and  his  co-workers^^"^^  developed  a  method  for  computing 
massively  separated  unsteady  incompressible  flow  fields.  This  method  was  based  on 
the  velocity-vorticity  formulation  of  the  Navier-Stokes  equations  which  consists  of 
the  vorticity  transport  equation  and  an  integral  equation  for  velocity.  The  results 
were  consistent  with  experimental  data. 

Cebeci,  et  have  reported  extensive  investigations  using  an  interactive 

approach  which  solves  inviscid  and  boundary-layer  equations  and  allows  them  to 
influence  each  other  in  an  iterative  manner.  This  method  was  developed  and 
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tested  for  steady  flows  and  was  used  in  a  quasi-steady  manner  to  examine  the 
evolution  of  the  flow  behavior  around  oscillatory  airfoils  operating  in  light  stall 
conditions.^^ 

Visbal  and  Shang^^"^^  presented  a  series  of  papers  describing  the  development  and 
application  of  a  compressible  Navier-Stokes  solver.  These  papers  address  many  of 
the  physical  variables  associated  with  dynamic  stall.  Rumsey  and  Anderson^® 
applied  an  upwind-biased,  implicit  approximate  factorization  algorithm  to  several 
unsteady  flows  on  dynamic  meshes.  They  used  the  thin-layer  form  of  the 
compressible  Navier-Stokes  equations  to  solve  both  laminar  and  turbulent  flows 
over  airfoils  pitching  about  the  quarter  chord.  Shida  and  Kuwahara^^  modeled  the 
time-accurate  static  stall  of  a  NACA  0012  airfoil  with  artificial  viscosity  that 
permitted  resolution  of  a  small  scale  structure.  Shida,  et  al?-^  further  modeled  the 
dynamic  stall  of  the  NACA  0012  airfoil  using  a  time-accurate  unsteady  Navier- 
Stokes  equation  solver  and  computed  the  flow  over  the  NACA  0012  airfoil 
oscillatory  in  pitch  at  M  =  0.3,  Re  =  4  x  10^ . 

Ono^^  simulated  the  dynamic  stall  process  on  a  two-dimensional  NACA0012  airfoil 
oscillating  in  pitch.  The  qualitative  agreement  with  experimental  data  was  fairly 
good,  but  quantitative  agreement  wasn't  satisfactory.  Currier  and  Fung^^  conducted 
a  numerical  study  to  assess  the  sensitivity  of  the  separating  boundary  layer  to  the 
transition  location.  They  found  that  the  bursting  of  the  separation  bubble  at  the 
airfoil  leading  edge  is  the  onset  mechanism  for  most  of  the  dynamic  stalls.  Jang  et 
al?-^  applied  the  implicit  approximate  factorization  solution  algorithm  of  Beam- 
Warming^^  to  the  computation  of  the  unsteady  boundary  layers  on  a  rapidly 
pitching  NACA0012  airfoil  and  found  good  agreement  with  the  unsteady  pressure 
measurements  of  Landon^^. 
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The  most  recent  computational  studies  have  been  presented  by  Shrewsbury  and 
Sankar^®,  Grohsmyer  et  al.^,  and  Ghia  et  al?^ 

1.2.2  Experimental  Approach 

Experimental  data  has  provided  the  basis  for  the  current  understanding  of  the 
dynamic  flowfield.  Most  of  the  experimental  data  for  unsteady  separated  flow  over 
airfoils  have  been  obtained  for  oscillatory  airfoils  undergoing  relatively  small 
sinusoidal  pitch  oscillations  (±1  -  ±10  degrees)  about  a  relatively  small  mean  angle  of 
attack  (0  - 15°)  as  typified  by  the  experiments  reported  by  McCroskey  and  Philippe, 
McAlister  and  Carr,^^  and  Martin,  et  al.P.  A  limited  amount  of  experimental  data 
have  been  obtained  for  airfoils  undergoing  constant  pitching  rate  motions  up  to 
moderate  angles  of  attack  of  at  least  30°.  These  works  include  the  study  of  Harper 
and  Flanigan,^^  who  obtained  force  balance  data  on  a  small  aircraft  model  pitching 
up  to  30°,  the  work  of  Ham  and  Garelick,^^  who  measured  surface  pressure  on  an 
airfoil  pitching  up  to  30°,  and  the  work  of  Francis,  et  al?^  who  measured  surface 
pressure  on  an  airfoil  pitching  up  to  60°.  None  of  the  above  mentioned  works 
contain  any  flow  visualization  data.  Deekens  and  Kuebler^^  obtained  flow 
visualization  data  from  a  NACA  0015  airfoil  and  observed  the  dynamic  leading-edge 
separation  phenomenon  for  several  low  Reynolds  numbers  (less  than  3  x  10^)  and 
non-dimensional  pitch  rates  up  to  0.26.  Daley^  obtained  leading-edge  dynamic  stall 
data  for  Reynolds  numbers  up  to  3  x  10^  and  non-dimensional  pitch  rates  up  to  0.06. 
Walker,  et  al.^^  obtained  flow  visualization  data  along  with  hot  wire  data  on  a 
NACA  0015  airfoil  undergoing  constant  pitch  rate  motions.  These  data  were 
obtained  for  a  Reynolds  number  on  the  order  of  4.5  x  10“*  and  non-dimensional  pitch 
rates  up  to  0.30. 

Compressibility  effects  have  been  addressed  experimentally  during  the  last  few 
years.  Results  of  Schlieren  studies  by  Chandrasekhara  and  Carr‘*°  on  an  oscillatory 
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airfoil  have  indicated  that  compressibility  effects  set  in  at  M  =  0.3.  Further  studies  by 
Chandrasekhara  and  Ahmed^^  using  LDV  have  shown  the  formation  of  a 
separation  bubble  near  the  leading  edge  prior  to  the  formation  of  a  dynamic  stall 
vortex.  Studies  using  the  PDI  technique  by  Carr,  et  have  confirmed  the 

presence  of  the  separation  bubble  and  have  shown  that  the  flow  gradients  are  slow 
to  develop  in  the  oscillatory  case  compared  to  the  steady-state  resulting  in  the  delay 
of  stall  known  as  dynamic  stall.  Ahmed  and  Chandrasekhara‘S^  carried  out  a 
detailed  study  of  the  reattachment  process  of  dynamic  stall  flow  over  an  oscillatory 
airfoil.  They  have  found  that  reattachment  progresses  through  a  separation  bubble, 
which  change  size  during  the  process  and  disappears  at  a  low  angle  of  attack. 

Chandrasekhara  and  Brydges'*^  documented  the  effects  of  increasing  amplitude  on 
an  oscillatory  airfoil  in  both  compressible  and  incompressible  flows  and  showed  that 
larger  amplitudes  resulted  in  vortex  retention  at  higher  angle  of  attack  for  a  given 
Mach  number  and  reduced  frequency. 

1.2.3  Three-Dimensionality  of  Flow 

The  effect  of  three-dimensionality  on  aircraft  dynamic  stall  is  significant  and  must 
be  included  at  the  onset  if  a  full  representation  of  dynamic  stall  on  an  aircraft  is  to  be 
attained.  The  effect  of  pitch  oscillation  on  a  delta  wing  was  studied  experimentally 
by  Gad-el-Hak  and  Ho.^^  They  found  significant  interaction  between  the  vortices 
shed  from  the  leading  edge  and  those  shed  during  the  dynamic  stall  process.  They 
also  studied  the  low  Reynolds  number,  and  time-dependent  flows  around  the  delta 
and  swept  wings'^^  and  found  that  on  the  rectangular  wing,  the  leading  edge 
separation  vortex  convects  downstream,  while  it  is  stationary  during  part  of  the 
angle  on  the  swept  wing.  On  the  delta  wing,  the  leading  edge  vortex  does  not 
convect,  rather  it  experiences  a  grow-decay  cycle. 
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Adler  and  Luttges'*®  obtained  flow  visualization  for  a  rectangular  wing  with  an 
aspect  ratio  of  two,  and  showed  that  flow  features  were  similar  to  two-dimensional 
stalls  appearing  within  a  chord  length  of  the  tip.  Carta studied  the  effect  of  sweep 
on  an  oscillatory  airfoil  and  found  that  sweep  effects  appeared  near  the  leading  edge, 
and  that  there  were  large  phase  shifts  in  the  lift  results  for  the  swept  and  unswept 
wing,  but  only  if  dynamic  stall  had  occurred  in  the  cycle.  Wagner^®  observed  that 
the  location  and  size  of  the  tip-vortex  changes  significantly  with  small  variation  in 
the  amplitude  of  oscillation. 

In  their  measurement  of  unsteady  pressure  distribution  on  a  pitching  rectangular 
wing,  Robinson  and  Wissler^^  observed  the  interaction  between  the  dynamic  stall 
vortex  and  the  tip  vortex  resulting  in  prolonging  the  passage  of  the  stall  vortex 
which  in  turn  enhances  the  value  of  sectional  lift  coefficient.  Ashworth,  et  al.^^ 
studied  three-dimensional  flow  field  about  a  forward  swept  NACA  0015  wing.  They 
found  that  strong  helical  tip  flow  vortices  dominated  most  of  the  observed  flow 
structures  near  the  wing  tip  across  all  test  conditions,  and  the  far  inboard  span 
locations  were  dominated  by  flows  related  to  the  leading  edge  vortex.  St.  Hilaire,  et 
al  53,54  examined  the  effect  of  sweep  on  an  oscillatory  wing  model.  They  found  that 
sweep  tends  to  delay  the  onset  of  dynamic  stall  and  slightly  reduces  the  magnitude 
of  the  hysteresis  loop.  Garta^^  found  that  near  the  leading  edge  the  sweep  effect  is 
significant,  and  there  are  phase  shifts  in  the  aerodynamic  forces  between  swept  and 
unswept  wings  when  dynamic  stall  has  occurred  in  the  cycle. 

Ashworth,  et  Luttges  and  M.C.  Robinson^®,  and  Adler  and  Luttges®^  have 

made  a  series  of  studies  on  three-dimensional  vortex  flows  created  by  sinusoidal 
oscillation  of  wings. 

Salari  and  Roache^®  investigated  the  influence  of  sweep  on  the  deep  dynamic  stall  of 
a  rapidly  pitching  swept  wing  using  numerical  simulations.  They  found  that  sweep 


tends  to  delay  the  onset  of  dynamic  stall  and  reduce  the  magnitude  of  unsteady 
aerodynamic  loads.  Chaderjian  and  Guruswamy^^  used  a  zonal  grid  approach  to 
simulate  unsteady  transonic  viscous  flow  about  a  three-dimensional  rectangular 
wing  with  an  oscillatory  angle  of  attack.  Computed  real  and  imaginary  pressure 
coefficients  compared  well  with  the  experimental  values. 

1.3  Phase  1  Study  and  Its  Merits 

In  this  project  we  developed  and  demonstrated  the  high-order  accurate  and  efficient 
pressure-based  algorithm  for  predicting  the  quantitative  features  of  separating 
unsteady  flows,  such  as  dynamic  stall  in  the  two-  and  three-dimensional  stationary 
and  maneuvering  bodies  of  aerospace  vehicles.  The  major  innovations  and 
contributions  to  the  state-of-the-art  of  the  Computational  Fluid  Dynamic  of  the 
pressure-based  method  are: 

1.  Development  of  high-order  TVD  scheme  for  the  pressure-based 
algorithm; 

2.  Newton's  iteration  technique  for  fast  convergence;  and 

3.  Novel  concept  of  hyperbolic  pressure  correction  equation  to  improve 
solver  efficiency. 

1.3.1  Why  Pressure-Based  Algorithm  ? 

Most  aerospace  vehicles  tend  to  operate  in  the  transonic  regime  where  the  flow  field 
is  primarily  subsonic  with  regions  of  supersonic  flows.  This  often  gives  rise  to 
complex  fluid  physics  such  as  steady  and  lime  dependent  vortical  flow,  shocks  and 
separations.  A  solution  algorithm  is  needed  which  is  equally  effective  for  a  high¬ 
speed  compressible  flow  regime  as  well  as  for  a  separated  and  stagnated  mildly 
compressible  regime. 


Current  state-of-the-art  of  CFD  technology  can  be  divided  into  two  groups: 
density-based  methods  and  pressure-based  methods.  In  the  density  based  approach, 
density  is  treated  as  a  transport  variable  in  the  continuity  equation.  Pressure  is 
derived  from  the  equation  of  state.  The  density-based  algorithms  have  been  widely 
used  for  compressible  and  external  flows.  The  merit  of  this  approach  lies  in  its 
ability  to  obtain  high  order  accuracy,  which  is  accomplished  by  applying  recently 
developed  high  resolution  schemes,  such  as  ENO^^'^®,  MUSCL^^  and 

PPM^o.  Indeed,  the  density-based  method  incorporates  the  non-linear  wave 
properties  into  the  numerical  solutions  in  the  form  of  Riemann  problems  and 
characteristic  equations.  This  leads  to  algorithms  which  are  robust  and  accurate  for 
high-speed  compressible  flows^^ 

The  accuracy  and  efficiency  of  the  density-based  methods,  however,  breaks  down  at 
low  Mach  numbers  and  for  recirculating  flows.  Here  acoustic  wave  speed  becomes 
very  high  relative  to  the  fluid  velocity,  and  CFL  restriction  requires  a  very  small 
time  step.  Whereas  in  the  incompressible  flow  limit,  density  is  constant  and  is 
independent  of  pressure,  so  that  the  pressure  field  in  the  momentum  equation 
cannot  be  extracted  from  the  density  field,  and  density-based  methods  fail.  To  date, 
the  methods  have  rarely  been  used  for  incompressible  or  low-speed  turbulent  and 
recirculating  flows.  Artificial  compressibility  has  to  be  introduced  in  the  continuity 
equation  to  use  the  density  based  approach  for  incompressible  flows. 

Pressure  based  methods,  on  the  other  hand,  are  effectively  characterized  by 
combining  the  continuity  and  momentum  equations  to  form  a  Poisson-like 
equation  for  pressure  or  pressure  correction.  Here  any  change  in  density  is  then 
considered  a  function  of  pressure  via  an  equation  of  state.  As  a  result,  it  can  handle 
both  compressible  and  incompressible  flows  with  equal  accuracy  and  efficiency.  This 
approach  has  been  very  successful  in  complex,  recirculating  and  turbulent  flows. 


10 


The  governing  equations  in  the  pressure-based  approach  are  usually  solved  in  a 
segregated  manner  (one  variable  at  a  time).  Instead  of  solving  block  matrix 
equations  in  a  factorized  form,  as  in  the  density  based  methods,  a  single  equation 
matrix  is  solved  on  the  entire  computational  domain.  For  elliptic  flow  problems,  as 
in  a  dynamic  stall  case,  it  is  very  efficient  and  requires  less  computational  storage. 
Recent  assessment  of  the  pressure-based  algorithms'^'  for  one-dimensional  fast 
transient  and  resonant  compressible  flows  with  shocks,  shows  a  very  promising 
future  for  the  method  in  high-speed  flows.  In  Figure  1-2  results  of  the  compressible 
flow  with  moving  shock  in  a  resonant  pipe  calculated  by  density  and  pressure  based 
methods  are  displayed.  The  selected  test  case  is  very  challenging  in  that  it  requires 
the  numerical  method  to  be  non-dissipative  and  non-dispersive,  be  able  to  capture 
shock  without  wiggles,  and  be  able  to  keep  shock  amplitude  for  a  long  time.  As  seen 
from  the  Figure  1-2,  the  proposed  pressure-based  method  is  as  accurate  as 
density-based  method. 

1.3.2  Merits  of  the  Present  Methodology 

There  are  several  merits  of  the  presently  proposed  methodology.  First,  the  present 
pressure-based  approach  has  the  advantage  of  being  a  unified  methodology.  It  can 
be  applied  to  a  wide  variety  situations:  unsteady  and  steady  flows;  low  speed 
subsonic,  transonic,  supersonic  and  hypersonic  flows;  perfect  and  real  gas;  viscous 
and  inviscid  flows;  single  and  multiple  spatial  dimensions;  simple  and  complex 
geometries;  internal  and  external  flows;  etc.  Secondly,  since  the  convective  terms 
in  the  governing  equations  are  modeled  by  a  high  order  TVD  scheme,  it  requires  no 
user-specified  dissipation  terms.  In  contrast,  most  of  the  existing  pressure-based 
codes  require  user-specified  dissipation,  which  can  result  in  lost  accuracy  of  the 
solution  not  only  near  discontinuities  but  also  across  the  computational  domain. 
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C-O  Third-Order  TVD  Scheme,  Pressure-Based  Method 


Figure  1-2.  Reflection  and  Propagation  of  Shock  Wave  in  One-Dimensional 
Tube  Calculated  by  Density-  and  Pressure-Based  Methods 
C-O  =  Chakravarthy-Osher  Third-Order  Scheme^^ 


Thirdly,  by  strongly  coupling  the  velocity  and  pressure  fields,  the  presently 
developed  Newton's  iteration  technique  greatly  improves  the  numerical 
convergence  rate  and  does  not  require  the  difficult  evaluation  of  Jacobian  matrix. 
Finally,  proposed  innovative  hyperbolic  form  of  the  pressure  correction  equation 
can  be  solved  more  efficiently  than  the  commonly  used  elliptic  equation  form. 

1.4  Technical  Objectives  and  Approach 

The  objective  of  the  proposed  project  is  to  develop  a  high  order  TVD  scheme  and 
efficient  pressure-based  algorithm  and  to  demonstrate  their  capability  for  a  model 
problem  of  stall  flow  on  static  and  dynamic  maneuvering  three-dimensional  wings 
in  subsonic  condition.  The  specific  objectives  of  Phase  I  effort  include: 

•  Develop  a  conservative  and  consistent  formulation  of  third-order  TVD 
scheme  applicable  for  conservative  and  primitive  variables; 

•  Develop  Newton's  iteration  technique  for  a  complete  set  of  primitive 
variables  (u,  v,  w,  p,  H,  k,  e,  ...)  solved  by  a  "whole-field"  rather  than 
Block-TDMA  type  or  Block-Gauss-Seidel  type  solvers.  Utilize  the 
maximum  available  information  from  Riemann  solutions  and  entropy 
condition; 

•  Perform  and  evaluate  the  effectiveness  of  hyperbolized  form  of  pressure 
correction  equation;  and 

•  Use  the  modified  code  to  investigate  the  flow  physics  of  dynamic  stall  on  a 
static  and  dynamic  maneuvering  three-dimensional  wing  in  subsonic 
flow. 

The  proposed  methodology  has  been  partially  evaluated  on  a  one-dimensional 
nonlinear  acoustic  problem  with  excellent  results^^.  The  challenge  for  this  project  is 
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to  extend  it  to  multidimensional  flows  and  retain  the  high  numerical  accuracy  for 
low  and  high  Reynolds  numbers,  and  for  sub-  and  supersonic  flows. 


The  development  of  the  above  new  methodologies  has  been  achieved  by  adapting 
and  modifying  an  existing  advanced  CFD  code,  REFLEQS.^^^^ 

The  REFLEQS  code  has  been  developed  by  CFDRC  personnel.  Its  main 
capabilities/methodologies  directly  related  to  the  present  work  are: 

1.  Solution  of  two-  and  three-dimensional  Navier-Stokes  equations  for 
compressible  and  incompressible  flows; 

2.  Cartesian,  polar,  and  non-orthogonal  body-fitted-coordinates  and  body- 
conforming  moving  grid; 

3.  Fully  implicit  and  strong  conservation  formulation;^^ 

4.  Central  differencing  with  damping  terms; 

5.  Second-order  time  accurate  formulation; 

6.  Four  turbulence  models:  Baldwin-Lomax  model;  Standard  k-e; 
multiple-scale  model  of  Chen;^  and  Low  Re  k-e  model  of  Chien; 

7.  Symmetric  and  periodic  whole  field  solver; 

8.  Pressure-based  solution  algorithms,  with  an  enhanced  variant  of 
SIMPLEC,  SIMPLE,  and  PISO;  and 

9.  User  friendly  pre-processor  and  graphical  p>ost-processor. 

Significant  emphasis  has  been  placed  on  the  systematic  and  quantitative  validation 
of  REFLEQS.  It  has  already  been  validated  for  over  thirty  benchmark  problems. 
References  73  through  75  describe  some  of  these  problems. 
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1.5 


The  next  section  (Section  2)  describes  the  basic  pressure-based  methodology  in 
current  REFLEQS  code.  Section  3  presents  the  proposed  improvements  including 
high-order  TVD  schemes,  Newton's  iteration,  and  hyperbolization  of  pressure 
correction  equation. 

Section  4  presents  the  validation  of  the  REFLEQS  code  on  2D  airfoils  for  inviscid 
supersonic  and  transonic  flows,  and  viscous  transonic  flow.  Several  dynamic  stall 
conditions  are  also  simulated.  Section  5  describes  the  simulations  and  discusses 
flow  physics  of  static  and  dynamic  stalls  on  three-dimensional  rectangular,  forward 
swept,  backward  swept ,  and  delta  wings.  Comparisons  are  made,  wherever  possible, 
with  experimental  data  and  visualization. 

Finally,  Section  6  presents  the  summary  of  the  present  study  and  recommendations 
for  further  investigations.  Preliminary  recommendations  for  the  SBIR  Phase  II 
study  are  also  outlined.  Further  selection  and  elaboration  of  these 
recommendations  will  be  presented  in  the  proposal  for  Phase  n  work. 


1.5 


2. 


PRESSURE-BASED  METHODOLOGY 


The  objective  of  the  present  study  is  to  develop  advanced  pressure-based 
methodology  to  study  the  physics  of  three-dimensional  dynamic  stall.  To 
accomplish  the  above  objective,  an  existing  advanced  CFD  code,  REFLEQS,  will  be 
modified  by  implementing  the  proposed  techniques.  This  section  will  briefly 
describe  the  basic  pressure-based  methodology  available  in  REFLEQS.  The  proposed 
more  advanced  techniques  and  their  uniqueness  will  be  presented  in  the  next 
section. 


2.1  Governing  Equations  and  Transformation 


The  flow  governing  equations  are  the  compressible,  Reynolds-averaged,  Navier- 
Stokes  equations.  They  can  be  written  in  a  Cartesian  tensor  form  as: 


dt 


+ 


dXj  3  dx^  '^1 


(2.1) 


(2.2) 


dt 


+  ^pu.W)  = 

dXj 


+  A 

till 

dt  dxj 

.  3xj  dxj 

^ax.  ax. 

dxi 


(2.3) 


where  p  is  density,  Uj  is  the  Cartesian  velocity  component,  p  is  the  pressure,  H  the 
total  enthalpy,  T  the  temperature,  and  t  the  time.  Xj  is  the  Cartesian  coordinate. 


=  X  ,  Xj  =  y  ,  X3  =  2 


(2.4) 
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|i  and  r  are  the  effective  dynamic  viscosity  and  effective  thermal  conductivity, 
respectively: 


H  =  +  (2.5) 

+  (2.6) 

The  subscript  1  denotes  laminar  quantity,  and  t  the  turbulent  quantity. 

For  most  engineering  applications,  including  the  present  problem,  Cartesian 
coordinates  are  rarely  adequate  in  describing  the  geometric  configuration.  Thus,  a 
generalized  coordinate  mapping  is  introduced  in  the  form 

^  /  y  =  y  /  2  =  2  ^  =  2  (2.7) 

or  simply 

=  ,  t  =  x  (2.8) 

where 

=  C  (2-9) 

The  purpose  of  introducing  temporal  variation  of  the  coordinates  is  to 
accommodate  body-conforming  moving  coordinates. 

Using  the  chain  rule  we  have. 


8  _  8  (2.10a) 

8t 
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(2.10b) 


A 

dxi  dxi  3^; 

By  some  manipulation.  Equations  2.1  to  2.3  can  be  written  in  a  strong  conservation 
form,  except  pressure  term: 


£ 

dt 


=  0 


a  ipui 

1  a  IpUjuA  Tdi’dp  , 

dt\  j 

a 

M  1 

d^^dui  d^^duj  2  dui  a.^^^ 

J  dXj  \ 

3*1 3^'  3  'dr  3*1/ 

(2.11) 


(2.12) 


d  |p//\  ^  9  IP^/^\  l^P  ,  ^  9^^  BT 

dr  J  ^  I  ’-^dt  d^’^ydxj  dxj^^‘i 


+ 


[ax.  3^'  8x^3^'  3  3^^^ 


(2.13) 


where  Uj  is  the  velocity  component  in  the  direction,  also  known  as  the 
contravariant  velocity  component. 


dt  dx^ 


Wi 


(2.14) 


X  has  been  replaced  by  t  in  the  above  equations.  J  is  the  coordinate  transformation 
Jacobian,  given  by 


(2.15) 


1^ 


_ 


idx2 


dx^  dx2 


a(xW) 


dx^  Idx^  dX2  BX2  9^:3  \  33C, 

a^'a^7  ar 


faX2  Bx^  Bxj  BX2 


13«’ 


a^’ai^ 


or 

j- = ^«  V?)  *  ^  -  yi"?) + (y«^i  -  J 

The  coordinate-transformation  matrices  are  defined  as 


^X  =  J  (yrjZ^  -  Zr,y^)  ,  Vx  =  J  '  y^f) , 
4=  /  {ZrjX^  -Xr,Z^),  r]y  =  }  {X^^  -  Z^f)  , 

^x  =  J{xr,y;-ynX;) ,  Vz=J{y^r^^d ' 

Cx=^J  [y^n  -  2^»j)  /  =  -  xt^x  •  yt^y  -  ztiz  / 

Cy  =  ;  (2|Z;,  -  Z|2^)  ,  m  =  -Ztr;^  -  ytVy  -  ZtHr  / 
Cz  =  J  {x^rj  -  y^n) '  Ct  =  "XtCx  -  ytCy  -  2tCz 


(2.16) 


(2.17) 


Unlike  the  density-based  method,  where  pressure  is  directly  expressed  in  terms  of 
density  and  kinetic  energy  by  the  equation  of  state,  the  pressure-based  method  keeps 
the  pressure  gradient  terms  in  the  momentum  equation  for  the  coupling  with 
velocity. 

2.2  Discretization  of  Governing  Equations 

The  discretized  finite-difference  forms  of  the  governing  differential  equations  are 
obtained  using  a  finite-volume  approach.  First,  the  solution  domain  is  divided  into 
a  finite  number  of  discrete  volumes  or  "cells,"  where  all  variables  are  stored  at  their 
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geometric  centers.  In  the  current  REFLEQS  code,  a  non-staggered  grid  system  is 
adapted,  and  a  typical  control  volume  is  shown  in  Figure  2-1. 


Equations  2.11  to  2.13  are  then  integrated  over  all  the  control  volumes  by  using  the 
Gauss's  theorem.  For  example,  the  integration  of  the  continuity  equation  of  2.11 
leads  to 


(P  Po^)p-(P  Hp  ,  r 
At  * 


+  G;,  -  G,  =  0 


(2.18) 


where  vol  is  the  volume  of  the  cell,  and  the  G's  represent  the  mass  flux  over  the 
control  volume  surfaces. 


/ 


/ 


'w 


C*  =  (^p4.  G,  =  (^pU3)_ 


(2.19) 


The  subscripts  e,  w,  n,  s,  h,  and  1  represent  the  values  at  east,  west,  north,  south, 
high,  and  low  faces.  The  superscript  "o"  represents  the  previous  time  level.  For 
clarity,  the  first-order  backward  Euler  differencing  scheme  is  shown.  For  second 
order  accuracy  in  time,  the  fluxes  of  G's  can  be  interpolated  between  the  present 
time  level  and  the  previous  time  level  "o".  The  detailed  derivation  will  be  omitted 
here. 
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Figure  2-1.  Curvilinear  Coordinates  and  Finite- Volume  Representation 
Integration  of  momentum  equation  2.12  results  in 

At  I  °  k 

I  i, 

^  ^  '  35‘1 
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where 


(2.21) 


dXj  dxj 


and 


S.  = 


dUj 

1  '  1 

\dx.  dx,.) 

2  l\ 

3^  3<^| 

3«^. 

1  ’  1 

\3xj  Sx,l 

ai' 

dUj 

1  '  1 

\dxj  Bx,l 

ai' 

.1  3 


ar1 

dUf 

dXil 

ar 

2 

a^ 

a<^ 

dUi 

\dXi 

dxi 

ar 

3 

wi  \ 

a^ 

a^ 

3m, 

\Bxj 

dXi 

a| 

I 


(2.22) 


6:; 


Similarly,  the  energy  equation  can  be  written  as; 


{p  vol  H)p  -  (p  vol  H)p 
_  + 


.2  .3 


^  ai‘j 


3 


K 


(2.23) 


-  Sfj 

In  the  above  equation,  the  heat  diffusion  terms  have  been  replaced  by  the  diffusion 
of  total  enthalpy.  is: 


22 


As  seen  from  above,  the  momentum  equation  and  the  energy  equation  are  in  the 
same  form.  The  key  issue  is  to  approximate  convective  terms  and  diffusion  terms 
from  cell  faces  (e,  w,  n,  s,  h  and  /)  to  cell  centers  (P,  E,  W,  N,  S,  H,  and  L).  We  will 
take  the  east  face  as  an  example.  The  convective  term  C  can  be  generally  written  as: 


Cg  —  ;  (p  —  Ui ,  H 


(2.25) 


The  following  approximation  can  be  made; 

=  +(l-p)d;^  (2.26) 


where  superscript  UP  stands  for  upwind  differencing  and  CN  for  central 
differencing. 


(2.27a) 


C“  =  fk-#p)  G.27b) 

As  a  result  Equation  (2.26)  becomes. 

As  for  the  diffusion  term,  D^: 

\  '  a« ;« 

= Ln  if. + ji2  il + „i3  il 

/  \  de  Bn  (2-29) 

_rA;^A|^/  fe-# 

/  r  A?'  *  A^^  A^^  I 


A  further  approximation  is  made  as  to  express  (p^^,  (Pf^g,  and  in  terms  of  the 
values  at  cell  center: 


te  =  li<pN+h^^p-^^NE) 

^se  =  4(^S  +  ^£  +  ^P  +  ^SE) 


(p^+  <Pp+  (pyi^) 


=  4  (^L  ■*■  <**£■*■  ^P  hp) 


(2.30) 


The  same  procedure  applies  to  the  convective  and  diffusive  terms  at  west,  north, 
south,  high,  and  low  faces.  With  these  representations,  a  finite-difference  analog  of 
the  governing  equations  can  be  expressed  as: 
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Ap  (l>p  =  Ap<l>p+  Ayf  0VV  +  h  \  h  + 

^SW^SW^  ^NW  'PnW  ■*■  ^SE  ^SE  ^NE  ^NE  ■*■ 

(2.31) 

^WL  ^WL  ^EL  ^EL  '*'  ^WH  ^WH  '*'  ^EH  ^EH  ■‘‘ 

^LS  ^LS  +  ^HS  ^HS  **■  \n  ^LN  ^HN  ^HN  '*'  ^ip 

or  simply 

Ap(t)p  =  X  A„y<l)„i,  +  X  ^cnb^cnb  +  (2-32) 

nb  cnb 

where 

0  =  Up  H 

'  nb=  E,  YJ,  N,  S,  H,  L 

cnb  =  SW,  NW,  SE,  NE,  WL,  EL,  WH,  EH,LS,  HS,  LN,  HN 

cnb  means  cross-neighbor.  A^^f^  vanishes  when  the  grids  are  orthogonal,  and  it  is 
small  compared  to  the  other  part  if  the  grid  non-orthogonality  is  not  severe. 

2.3  Pressure-Velocity  Coupling 

The  discretized  finite  difference  equation  (2.32)  does  not  apply  to  continuity 
equation.  Unlike  the  density-based  method  where  density  is  treated  as  a  transport 
variable  in  the  continuity  equation,  the  pressure-based  method  takes  pressure  as  the 
dependent  variable  in  the  continuity  equation.  It  requires  a  proper  coupling 
between  velocity  and  pressure.  Presently,  the  coupling  of  pressure  and  velocity  is 
achieved  via  the  well  known  SIMPLE  algorithm^®  and  its  variants  SIMPLEC^^  and 
PISO.®0 

As  seen  from  Equation  2.18,  the  contra  variant  velocity  component  is  needed  at  a 
control  volume  face.  For  example,  at  the  east  face. 
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[a«’| 

U'l ' 

(2.33) 


Since  only  the  values  L/^  at  cell  center  P  and  E  are  known,  an  interpolation  is 
required.  The  first  choice  may  be  the  linear  interpolation.  This,  however,  will  give 
rise  oscillations  in  the  final  resolution,  as  explained  by  Patankar^®.  Early  practice  was 
to  use  a  staggered  grid  arrangement  to  avoid  any  possible  oscillation. 

In  the  staggered  grid  arrangement,  velocity  components  at  cell  faces  are  solved  from 
momentum  equations  and  hence  no  interpolation  for  the  velocity  components  is 
necessary.  In  the  current  REFLEQS  code,  a  non-staggered  grid  system  is  used,  and  to 
avoid  oscillations,  a  special  interpolation  practice  is  employed,  as  suggested  by  Rhie 
and  Chow®^  and  by  Peric®^  . 

The  idea  is  based  on  the  supposition  that  if  the  velocity  values  at  the  cell  faces  were 
obtained  by  solving  the  momentum  equations  at  the  cell  faces,  these  momentum 
equations  would  contain  a  pressure  gradient  which  could  be  evaluated  as  the 
difference  between  neighboring  pressure  locations.  To  illustrate  this  procedure,  we 
write  the  momentum  equation  as: 


(»i)p  = 


^nb^nb  5-r  ^cnb^cnb  ^  ^Uj 


W-—] 

-P 

L  J 

(2.34) 


where 


Dii  = 


ApJ  Bxi 


(2.35) 
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In  comparison  to  Equation  2.32,  the  pressure  gradient  terms  have  been  taken  out  of 
and  shown  explicitly.  To  evaluate  (at  east  face)  the  terms  on  the  right-hand 
side  of  Equation  2.34  are  selectively  interpolated  for  "e"  location.  Thus, 


\  _  ^  ^nb^nb  ^  ^cnb^cnb  f-p 


(2.36) 


where  the  overbar  denotes  linear  interpolation.  The  cell  face  velocities  are  thus 
made  dependent  on  the  pressure  at  two  neighboring  nodes,  as  is  the  case  in  the  true 
staggered  arrangement.  For  uniform  grid,  with  the  use  of  Equation  2.34,  we  can 
write  the  above  expression  as: 


ai'lr  m'. 


^D.y 


(2.37) 


It  is  equivalent  to  linear  interpolation  with  added  third  order  pressure  damping. 
This  pressure  damping  term  is  necessary  for  incompressible  flows,  and  it  may  cause 
accuracy  problems  for  source  dominant  (such  as  rotational  and  gravitational)  flows, 
as  well  as  compressible  flows.  Therefore,  a  coefficient  is  multiplied  to  this  term: 


(«.)e  =  ^[(«f)E  +  («.)p]  +  ^ij - [(C>.y)E  +  (Dij)p]  +  Dij -^1 

\\  If  L  Ip 


(2.38) 


For  incompressible  flow  =  1.0,  and  for  compressible  flows  =  0.1.  With  [u-\ 
calculated  as  the  above,  the  contravariant  component  can  be  found  from  Equation 
(2.33). 
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From  now  on  the  SIMPLE  strategy  can  be  used  in  its  standard  form,  developed  for 
the  staggered  arrangement  grid  system.  The  mass  imbalance  resulting  from 
continuity  equation  2.18  is: 


Gj  +  G;,-Gj  = 


(2.39) 


where  superscript  *  denotes  values  employed  in  and  resulting  from  the  momentum 
equations. 


To  enforce  the  continuity  equation,  a  mass  flux  correction  G'  is  introduced,  which  in 
turn  relates  velocity  correction  u,-,  which  is  fxurther  related  to  the  pressure  correction 
p'hy: 

=  <2.40) 

«  V  I 

The  continuity  equation  then  reads 

^£^+g;-g„+g.-g;+g*-g;+s„=o 


For  Gg  we  have 


(2.42) 


From  equation  of  state  and  Equation  2.33,  we  have 


•  dp  . 


(2.43) 


(3«'l  („ 

■)  -H’] 

(2.44) 
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This  leads  to 


jdp 

I 

idp 

II ”1 

At 


3^,'  3|'/J 


2  3 

A^  A^ 

dp 

J 

[3p 

U,P'-P^D-^ 

3«' 


u; 


(2.45) 


l^V|ggu  p  .pgi!^..jgl 
^r]ap  '’a., 


=  -s 


m 


Equation  2.45  is  the  same  form  as  Equation  2.20,  and  hence  it  leads  to  an  equation  for 
the  pressure  correction  which  has  the  same  form  as  Equation  2.32. 


For  the  solution  domain  as  a  whole  this  results  in  a  system  of  N  equations  with  N 
unknowns,  where  N  is  the  number  of  control  volumes.  Several  efficient  iterative 
solvers  have  been  employed  to  solve  the  system  equations,  for  example,  modified 
strongly  implicit  procedure  (SIP)  of  Stone,®^  based  on  an  incomplete  LU  factorization 
of  the  coefficient  matrix. 


2.4  SflUttioiL.  Algorithm 

Based  on  the  proceeding  discretization,  the  solution  algorithm  for  the  transient 
flows  in  body-fitted  moving  coordinates  can  be  summarized  as  follows: 

1.  Specify  initial  grid  and  values  of  dependent  variables  (initial  conditions); 

2.  Determine  grid  velocity  from  the  new  grid  position  after  the  time 

has  advanced  by  At; 

3.  Calculate  the  coefficients  of  the  momentum  equation  and  solve  to  obtain  a 
new  velocity  field; 

4.  Calculate  new  values  of  mass  fluxes  through  the  cell  faces  using  Equations 
2.33  and  2.19,  and  determine  the  mass  imbalance  using  Equation  2.39; 
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5.  Calculate  the  coefficients  of  the  pressure-correction  equation  and  solve  to 
obtain  the  pressure-correction  field; 

6.  Correct  the  mass  fluxes,  nodal  velocity  components  and  pressure  by  the 
calculated  pressure; 

7.  Calculate  the  coefficients  of  other  scalar  equations  which  may  be  coupled 
with  the  momentum  equation  (e.g.  total  enthalpy,  turbulent  kinetic 
energy  and  its  dissipation  rate,etc.)  and  update  the  fluid  properties 
(density,  viscosity)  if  necessary; 

8.  Return  to  step  3  and  repeat  until  the  sum  of  the  absolute  residuals  in  the 
momentum  and  continuity  equations  has  failed  by  a  specified  order  of 
magnitude;  and 

9  Advance  the  time  by  another  increment  At  and  return  to  step  2;  repeat 
until  the  prescribed  number  of  time  steps  is  completed  or  the  prescribed 
time  has  been  reached. 

The  number  of  iterations  per  time  step  (steps  3-7)  depends  on  the  size  of  the  time 
increment  At;  for  smaller  At  fewer  iterations  are  needed  to  reach  the  solution  at  the 
new  time  level. 

2.5  Turbulence  Modeling 

When  the  Reynolds  number  is  high,  the  flow  around  an  airfoil  or  wing  is  turbulent. 
Due  to  the  finite  resolution  of  computational  grids,  a  computer  simulation  will  not 
be  able  to  capture  the  small  scale  vortices.  This  requires  the  modeling  of  the 
turbulence  effect.  In  the  current  REFLEQS  code,  there  are  four  turbulence  models 
available:  Baldwin-Lomax  model;  Standard  k-e;  Multiple  Scale  k-e,  and  Low 
Reynolds  number  k-e,  proposed  by  K.Y.  Chien. 
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Previous  studies  on  the  effect  of  turbulence  models  on  dynamic  stall  have  been 
conducted  by  Wu  et  al^.  They  compared,  for  a  variety  of  flows,  three  eddy  viscosity 
models:  the  algebraic  Baldwin-Lomax  model,  the  one  equation  Johnson-King 
model,  and  the  two  equation  Gorski  k-e  model.  They  found  that  the  Baldwin- 
Lomax  model®^  is  as  reliable  and  accurate  as  any  other  model  for  massively 
separated  flows,  but  that  all  models  are  to  date  inaccurate  in  the  same  regime. 
Furthermore,  the  Baldwin-Lomax  model  is  simple  to  employ  and  cost  efficient.  In 
light  of  Wu  et  al.'s  conclusions®'^,  we  will  simulate  the  effect  of  turbulence  by  the 
Baldwin-Lomax  two  layer  eddy  viscosity  model®^. 

2.6  Grid  Generation 

The  time-dependent  coordinate  transformation  (i.e.  moving  grid)  required  for  the 
present  flow  simulation  was  implemented  using  a  "rigid"  grid  attached  to  the  airfoil 
or  the  wing.  This  approach  eliminates  the  need  for  multiple  grid  generation,  and 
only  one  grid  is  required.  In  the  present  study,  an  O-H  grid  topology  is  employed, 
and  a  typical  three-dimensional  grid  around  a  wing  is  shown  in  Figure  2-2.  First  an 
O-grid  is  generated  using  algebraic  method  around  an  airfoil.  An  elliptic  solver 
based  on  Thompson  et  al.'s  methods®^  is  then  applied  to  improve  grid 
orthogonality,  if  required.  An  algebraic  method  is  then  used  to  generate  three- 
dimensional  grid  from  2-D  O-grid. 
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2.7 


The  present  procedure  requires  boundary  conditions  to  be  set  on:  a)  solid  boundary; 
b)  symmetric  boundary;  c)  far  field  boimdary;  and  d)  periodic  boimdary. 

2.7.1  Solid  Boimdary 

At  the  solid  boundary  the  "no-slip"  condition  requires  the  fluid  velocity  to  be  the 
same  as  that  of  the  solid,  and  the  motion  of  the  solid  is  a  known  function  of  time. 
Also  adiabatic  flow  conditions  have  been  applied  on  the  solid  surface.  The  pressure 
at  the  solid  surface  may  be  evaluated  in  various  ways,  such  as  solving  the  normal 
momentum  equation  at  the  surface.  In  the  present  study,  a  two  point  extrapolation 
was  iised.  For  example,  if  the  east  face  is  a  solid  boundary  and  uniform  grid  is  used: 

Pe  =  \Pp-\P\N  (2.46) 

The  above  form  was  found  to  give  virtually  the  same  results  as  those  obtained  by 
solving  the  normal  momentum  equations  at  the  solid  surface. 

2.7.2  Symmetric  Boundary 

This  boundary  applies  to  a  three-dimensional  wing  root  where  the  viscous  effect  is 
neglected.  It  also  applies  to  the  wall  for  invisdd  calculation.  At  this  boundary,  the 
normal  contravariant  component  is  set  to  zero,  and  other  components  are 
extrapolated  from  the  interior.  The  pressure  is  found  by  Equation  2.46. 
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2.7  J  Far  Field  Botindary 


At  the  far  field  boundaries,  except  in  the  downstream  boimdary,  the  flowfield  is 
assumed  to  be  prescribed  and  undisturbed.  At  the  downstream  boundary,  the 
velocity  field  and  total  enthalpy,  as  well  as  pressure,  are  extrapolated  from  the 
interior. 

2.7.4  Periodic  Boundary 

On  the  O-grid  cut,  spatial  periodicity  is  imposed.  This  leads  to  a  system  of  periodic 
matrix  equations. 
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3. 


PROPOSED  NEW  METHODOLOGIES  AND  IMPROVEMENTS 


The  pressure-based  numerical  solution  procedure  for  the  dynamic  stall  around  an 
airfoil  or  a  wing  has  been  described  in  the  previous  section.  In  comparison,  most  of 
the  previous  numerical  studies  relied  on  the  density  based  method.  The  typical 
representatives  of  density-based  codes  are  Pulliam-Steger's  ARC2D  and  ARC3D  flow 
solver®^"*®.  These  codes  solve  either  the  Euler  or  the  simplified  Reynolds-averaged 
Navier-Stokes  equations  with  the  standard  thin-layer  approximation.  Two 
numerical  algorithms  are  generally  supported  by  the  ARC  3D  flow  solver,  including 
ADI  algorithm  that  solves  block-tri-diagonal  matrices  along 
each  coordinate  direction  due  to  Beam  and  Warming^®  and  a  diagonalized  ADI 
algorithm  that  solves  scalar  pentadiagonal  matrices  along  each  coordinate  direction 
due  to  Pvilliam  and  Chaussee®^. 

In  comparison  to  the  density-based  method,  which  is  generally  limited  to 
compressible  flows,  the  pressure-based  method  can  solve  all  speed  flows  ranging 
from  incompressible  to  high  Mach  number  compressible  flows.  Instead  of  solving 
systems  of  equations,  the  pressure-based  method  solves  one  variable  at  a  time,  in  a 
segregated  approach,  which  requires  less  computer  storage.  Segregated  solution 
approach,  however,  may  lead  to  numerical  difficulties  when  strong  coupling  among 
all  of  the  variables  exists  {e.g.  high  Reynolds  number  flows).  The  presently  proposed 
Newton's  iteration  technique  aims  to  resolves  this  problem. 

High  resolution  TVD  schemes  have  been  well  proven  for  density-based  method, 
and  are  still  to  be  explored  for  pressure-based  method.  This  section  describes  three 
new  methodologies: 
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1.  High-Order  TVD  Schemes,  improving  the  accuracy; 

2.  Newton's  Iteration  Technique,  accelerating  numerical  convergence  rate; 
and 

3.  H5q)erbolic  Pressure  Correction,  for  pressure  solver  efficiency. 

3.1  High-Order  TVD  Schemes 

Total  Variation  Diminishing  (TVD)  concept  was  first  introduced  by  Harten^^.  An 
excellent  consolidation  of  various  TVD  methods  has  been  given  by  Yee^°.  These 
high  order  TVD  schemes  have  shown  excellent  accuracy  in  solving  Euler  equations 
for  compressible  inviscid  flows.  They  have  also  been  applied  in  a  straightforward 
manner  for  the  evaluation  of  the  convective  terms  of  the  Navier-Stokes  equations 
within  the  framework  of  density-based  algorithms.^'*®®  The  Research  Group  at 
CFDRC  has  initiated  a  study  on  high  order  TVD  schemes  in  the  pressure-based 
algorithm,^^  and  the  results  obtained  for  quasi-one-dimensional  nozzle  (Figure  3-1), 
ID  Riemann  problem  (Figure  3-2),  and  resonant  pipe  problem  (Figure  1-2)  are  very 
promising. 

The  fundamental  success  of  the  TVD  schemes  lies  in  the  accurate,  conservative  and 
consistent  evaluation  of  the  convective  flux  (Equation  2.26).  The  key  issue  is  how  to 
evaluate  the  "damping"  parameter  p  of  Equation  (2.26),  so  that  the  solution  is  high- 
order  accurate  and  oscillation  free.  In  TVD  schemes,  the  damping  p  is  adaptive 
based  on  the  profile  of  the  local  transport  variable.  By  manipulating  Chakravarthy 
and  Osher's  TVD  scheme^^,  we  can  find  the  east  face  flux  as; 

where 

p  =  minimod  minimod  (1 ,  cor) 


(3.1) 

(3.2) 
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with  6>  as  a  compressibility  parameter. 


(3.3) 


and  r  is  the  flux  difference  ratio: 


0)  = 


3-^ 


1-0 


r 


(c<=*'-c“’U 


,  and 


<T=  sign  {U^) 


(3.4) 


Figure  3-1.  Mach  Number  Calculated  Along 
a  Converging-Diverging  Nozzle  Length  by  Upwind, 
Central,  and  TVD  Schemes  With  Pressure-Based  Method 
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Figure  3-2.  1-D  Shock  Tube  Problem  by  TVD  Scheme  With  Pressure-Based  Method 

Based  on  the  value  of  0,  one  can  recover  different  schemes 

<j>  =  1/3  Third-order  scheme 

=  - 1  Fully-upwind 

=  0  Fromm's 

=  1/2  Low  truncation  error  second-order  (3.5) 

=  1  Central 

=  - 1/3  No  Name 

The  minimod  function  is  defined  as: 


minimod  (a,b)  =  sign  {a)  max  (o,  min  [w, 


(3.6) 


The  mathematical  meaning  of  this  limiting  procedure  is  to  pick  up  the  smallest 
value  between  a  and  b,  if  both  have  the  same  sign,  or  to  pick  up  zero  if  both  have 
different  signs. 


3.2  Newton^s  Iteration  Technique 

It  is  well  known  that  Newton's  iterative  method  gives  a  superior  convergence  rate. 
However,  Newton's  method  is  frequently  used  to  solve  small  systems  of  nonlinear 
algebraic  equations,  and  it  is  less  often  used  for  the  large  systems  of  nonlinear 
equations,  such  as  those  generated  by  the  discretization  of  Navier-Stokes  equations 
for  fluid  dynamics.  This  is  because  Newton's  method  requires  EXACT  evaluation 
and  inversion  of  the  Jacobian  matrix,  which  is  air  ost  impractical  due  to  limited 
computer  memory  and  speed.  Many  of  the  approaches  for  solving  Navier-Stokes 
equations  are  only  an  approximation  to  the  pure  Newton's  method.  For  example 
the  popular  Alternating  Direction  Implicit  (ADI)  method  in  density-based 
algorithms  uses  an  approximate  factorization  of  the  Jacobian  matrix  to  reduce 
numerical  operations  as  well  as  memory  storage  requirements.  In  addition,  it  is 
quite  common  to  make  approximations  for  elements  of  the  Jacobian  corresponding 
to  algebraically  complicated  terms  or  terms  that  increase  the  bandwidth  of  the 
factorized  Jacobian  matrix,  e.g.  terms  that  arise  from  the  turbulence  model. 

The  proposed  Newton's  method  requires  less  storage  while  making  no 
approximations  to  the  Jacobian  matrix.  To  illustrate  this  method,  we  use  Equation 
2.32: 


*  ^nb^nb  '  ^cnb^cnb  ~  (3.7) 

where  0  can  be  velocity,  total  enthalpy  or  pressure  correction.  The  nonlinearity  of 
the  above  equation  comes  from  the  link  coefficient  A's,  due  to 
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a.  velocity  quadratic  convective  features; 

b.  variable  density  in  convection  term;  and 
c  nonlinear  turbulent  diffusion. 


Let  subscript  "o"  denote  the  last  time  step  value.  The  standard  linearization  in  the 
pressure-based  method  is: 


(3.8) 


(f)  =(p~(j> 


The  convergence  rate  of  the  above  method  is  first-order  in  a  single  time  step.  We 
propose  a  Newton's  iteration  method  as: 


Ap(f>p  ~  ^  4r  Ap<Pp  +  A 

~  Ap<t>p  +  £  ^cnb^cnb 


(3.9) 


where  '*  means  the  last  iteration  value.  The  iteration  loop  is  imposed  for  velocity 
and  pressure  equations,  so  that  it  will  not  only  accelerate  the  convergence,  but  will 
also  strongly  couple  all  the  variables.  The  method  is  efficient  for  steady-state  as  well 
as  for  transient  problems. 


For  the  turbulence  model,  e.g.  k-e  model,  the  Newton  linearization  (Equation  3.9) 
also  applies  to  the  k  and  e  equations.  This  will  couple  the  whole  system  equation 
and  will  result  in  fast  convergence. 


3.3  Hyperbolic  Pressure  Correction 


It  has  been  recognized  that  the  pressure  correction  equation  (p'  equation)  in  the 
pressure-based  method  has  difficulties  to  converge,  specially  for  convection 
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dominated,  high  Reynolds  number,  incompressible  and  compressible  flows.  Even 
though  various  solvers,  such  as  block  correction,  conjugate  gradient,  coupled  whole 
field  solver,  multi-grid  method,  ...  have  been  tested,  success  is  still  limited.  Our 
experience  indicates  that  the  residual  error  of  momentum  equations  can  be  reduced 
to  machine  round-off  in  just  a  few  iteration  sweeps  of  linear  equation  solver. 
However,  the  iterative  convergence  rate  of  pressure  correction  equation  is  much 
slower  than  that  of  momentiun  equations.  It  is  particularly  difficult  to  converge  the 
p'-equation  for  fine  grids  with  large  cell  aspect  ratios  and  for  high  Reynolds  number 
flows.  We  believe  that  the  difference  between  momentum  equations  and  pressure 
correction  equation  in  the  finite  difference  form  lies  at  link  coefficients  (or 
coefficient  matrix).  For  the  momentum  equation,  when  convection  is  important, 
the  link  coefficients  are  hyperbolic.  In  the  sense,  for  one-dimensional  equation  if 
velocity  is  positive,  the  simple  upwind  scheme  gives: 


where 


D  =  diffusion  flux,  C  =  convective  flux 


(3.10) 


Whereas  for  the  pressure  correction  equation,  the  link  coefficients  are  always  elliptic 
regardless  the  importance  of  convection.  By  being  elliptic,  it  means  that 


A^-Dg,  (3.11) 

In  a  multi-dimensional  iterative  solver,  when  the  link  coefficients  are  hyperbolic, 
the  residual  (or  signal)  can  propagate  from  one  end  to  another  end  in  one  iteration 
sweep.  Whereas  for  elliptic  coefficients,  the  propagation  of  the  residual  is  by 
"diffusion,"  which  is  much  slower,  and  needs  a  large  number  of  iteration  sweeps. 
The  "elliptic"  nature  of  pressure  correction  equation  comes  from  the  general 
assumption:  velocity  variation  is  proportional  to  pressure  gradient: 
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This  relation  is  Darcy's  law,  and  it  holds  for  very  viscous  flow.  Actually,  by  using 
this  relation,  one  can  obtain  the  lubrication  equation  (Reynolds  Equation)  or  the 
equation  for  porous  medium.  On  the  other  hand,  for  invisdd  flow,  Bernoulli 
equation  says: 


^  pu^  +p  =  const  (3.13) 

or 

u'-p'  (3.14) 

if  we  substitute  this  relation  int ..  continuity  equation  for  incompressible  flow,  we 
have: 

Vm’  =  -Vu*  (3.15) 

which  is  hyperbolic  for  pressure  correction.  However,  if  we  substitute  Equation  3.12 
into  Equation  3.15  it  becomes  an  elliptic  equation.  Here  we  propose  an  innovative 
relation  for  pressure-velocity  relation  of  the  form: 

u  '^ap'  +  i>Vp’  (3.16) 

a,  in  the  above  equation,  is  a  coefficient  representing  hyperbolic  contribution,  and, 
b,  is  that  representing  elliptic  contribution.  To  determine  the  constants  a  and  b,  let 
us  consider  one-dimensional  inviscid  flow  of  the  form: 

8m  a(a2)_  Vp  (3.17) 

8r  dx  '  p 

After  discretization,  we  have: 


(3.18) 


where 


Sp  =  ^;  Vui  =  ui-u„b 


At 


(3.19) 


The  incremental  form  of  the  above  equation  is: 


(3.20) 


From  the  above  equation,  we  expect  the  following  asymptotic  expression: 


when^p^O, 


when  ->  0 ,  u,-  ~ 


VOl  Sr 


p(s^  + 


vp; 


(3.21) 


one  of  the  approximations  for  a  and  b  leads  to: 


U:  =  - 


vol  Anh 


P(sp  +  ^nj  *  p{4  +  A^i,) 


vol  Sr 


■Vpi 


(3.22) 


Now  the  pressure  correction  equation  should  contain  both  hyperbolic  and  elliptic 
coefficients.  One  should  be  able  to  reduce  the  residual  error  of  the  pressure 
correction  equation  easily  by  standard  solver,  without  using  multi-grid  or  block 
correction. 
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4.  TWO-DIMENSIONAL  STATIC  AND  DYNAMIC  STALLS 


The  Navier-Stokes  solver,  and  its  improvements,  described  in  the  previous  sections 
has  been  applied  to  a  number  of  test  studies,  for  the  purpose  of  validation  and 
studying  the  physics  of  unsteady  separation  flows.  In  all  cases,  the  computed  results 
have  been  compared  with  available  numerical  and  experimental  data.  This  section 
considers  the  following  test  problems: 

1.  Steady,  inviscid  transonic  and  supersonic  flows  past  a  NACA  0012  airfoil; 

2.  Steady  viscous  transonic  flow  over  a  NACA  0012  and  an  ARE  2822  airfoil; 

3.  NACA  0015  airfoil  undergoing  constant  rate  pitching  motions; 

4.  NACA  0012  airfoil  undergoing  oscillatory  pitching  motions;  and 

5.  Parametric  effect  on  dynamic  stall  on  the  airfoils. 

Each  of  the  above  cases  will  be  described  and  the  results  will  be  discussed  in  the 
following  sub-sections. 

4.1  Steady.  Inviscid  Transonic  and  Supersonic  Flows  Past  a  NACA  0012  Airfoil 
The  following  calculations  were  performed  for  the  NACA  0012  airfoil. 

a.  M^  =  1.2,  a  =  7.0° 

b.  =  1.8,  a  =  7.0° 

c.  M^  =  0.8,  0  =  1.25° 

These  cases  have  been  chosen  because  of  the  existence  of  accurate  numerical  data, 
especially  from  density-based  TVD  code,  available  for  detailed  comparison.  These 
cases  are  also  used  for  proving  the  superiority  of  proposed  TVD  methodology  and 
Newton's  iteration  technique. 


A  100  X  50  O-giM  is  used.  Figure  4-1  shows  the  grid  near  the  airfoil  surface.  The 
outer  boundary  is  placed  at  12c  distance  (c  is  the  chord  length). 


Figure  4-1.  Local  View  of  a  100  x  50  O-Grid  Around  a  NACA  0012 
Airfoil  for  Inviscid  Computations 

Figures  4-2  and  4-3  show  the  Mach  number  and  pressure  contours  around  the 
NACA  0012  at  a  =  7.0®,  =  1.2,  and  =  1.8  obtained  with  the  present  third-order 

TVD  scheme.  For  comparison  purposes,  those  from  Yee's  density-based  TVD,^  and 
the  widely  distributed  computer  code  ARC2D,  version  150^^  are  also  displayed. 
Their  results  are  obtained  using  a  163  x  49  C-grid.  It  is  seen  that  the  present  TVD 
scheme  gives  a  well-ordered  flow  structure  and  captures  the  shocks  with  a  good 
resolution.  It  is  very  competitive  with  the  density-based  TVD  scheme,  and  it 
actually  performs  better  at  the  trailing  edge  of  the  airfoil.  The  ARC2D  code,  on  the 
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Figure  4-2.  Comparison  of  the  Present  Pressure-Based  TVD,  Yee's  Density-Based  TVD 
and  ARC2D  for  Supersonic  Flow  Over  a  N ACAOOl  2  Airfoil  at  =1.2  and  a  =  7.0® 
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other  hand,  did  rather  poorly.  The  version  of  ARC2D^^  is  based  on  the  Beam- 
Warming  ADI  algorithm,^^  and  uses  a  mixture  of  second-  and  fourth-order 
numerical  dissipation  terms.  Due  to  the  adaptive  damping  property  inherent  in 
TVD  schemes,  one  can  capture  a  shock  in  one  to  two  grid  points  without  the 
associated  spurious  oscillations  near  shock  waves  as  opposed  to  three  to  four  by 
ARC2D. 

The  above  argviment  is  proved  in  Figure  4-4  which  shows  the  pressure  distribution 
around  the  NACA  0012  airfoil  when  =  0.8  and  a  =  1.25®  calculated  by  the  present 
TVD  scheme.  From  the  pressure  coefficient,  one  can  clearly  see  that  there  is  only 
one  transonic  point.  The  quality  of  the  result  is  remarkable. 

To  prove  the  superior  convergence  rate  of  the  proposed  Newton's  iteration  method, 
along  with  hyperbolic  pressure  correction.  Figure  4-5  shows  the  comparison  between 
the  regular  iteration  method  by  SIMPLEC  algorithm  and  the  present  Newton's 
iteration  method.  Undoubtedly,  residual  errors  drop  dramatically  at  a  constant  rate 
by  the  present  Newton's  iteration,  whereas  traditional  iteration  method  deteriorates 
and  slows  down  with  the  number  of  iterations.  Most  importantly,  for  some 
problems  where  shock  is  strong,  regular  iterative  method  cannot  converge,  while 
the  Newton  iteration  can  bring  the  residual  error  down  by  several  orders  of 
magnitude  without  any  difficulties! 

4.2  Steady.  Viscous  Transonic  Flow  Over  a  NACA  0012  and  an  ARE  2822  Airfoil 

A  Viscous  Transonic  Airfoil  Workshop  was  held  during  AIAA  25th  Aerospace 
Sciences  Meeting  in  1985.  In  this  workshop  a  series  of  two-dimensional  airfoil 
Navier-Stokes  computations  were  presented  to  establish  the  capability  of  various 
methods  for  computing  viscous  flow  fields  for  a  range  of  conditions  and  geometries. 
A  compendium  of  the  results  is  given  by  Holst.^^ 


48 


Iteration  Number 


Figure  4-5.  Convergence  History  of  Regular  Iteration  Method  and  Newton's 
Iteration  Method  with  Hyperbolic  Pressure  Correction  for  Invisdd  Flow 
Over  a  NACA  0012.  Grid  100  x  50,  M^  =  0.8,  and  a  =  1.25“ 


50 


The  purpose  of  the  present  simulations  is  to  evaluate  the  performance  of  the 
modified  REFLEQS  code  in  the  numerical  simulation  of  complex  turbulent  flows. 
The  test  cases  have  included  a  variety  of  different  situations  ranging  from  attached 
subcritical  to  transonic  flows  with  both  shock-induced  and  angle-of-attack-induced 
separations. 

The  Baldwin-Lomax  algebraic  eddy  viscosity  turbulence  model®®  has  been  employed 
in  the  present  study  for  the  analysis  of  transonic  airfoil  flow.  The  model  was  chosen 
for  its  computational  speed  and  simplicity.  In  addition,  it  has  been  shown  to  yield 
accurate  results  for  many  steady  flow  airfoil  computations.®®'®®  The  solutions  are 
assumed  to  have  reached  convergence  when  the  residual  of  all  equations  has 
dropped  at  least  five  orders  of  magnitude. 

4.2.1  NACA  0012  AirfoU 

Three  computations  were  made  for  the  NACA  0012  airfoil  corresponding  to 
experimental  results  of  Harris.^®®  All  three  were  computed  at  a  Reynolds  number  of 
9  million.  The  grid  used  for  the  present  NACA  0012  airfoil  is  a  200  x  63  O-mesh 
with  outer  boundary  extending  up  to  12c  (see  Figure  4-6  for  local  view  of  the  grid). 

Figure  4-7a  shows  the  computed  surface  pressure  coefficients  in  comparison  with 
experiments  at  =  0.7,  and  a  corrected  angle-of-attack  of  1.49®.  Agreement  is 
excellent.  For  this  case,  the  flow  is  attached  and  just  slightly  supersonic  near  the 
leading  edge  upper  surface.  From  density  contours  of  Figure  4-7b,  one  can  identify 
the  development  of  viscous  boundary  layer  on  the  airfoil  surface.  This  comparison 
indicates  satisfactory  performance  of  the  present  code  for  attached,  weakly  transonic 
flow. 
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At  a  higher  ai\gle-of-attack  of  8.34°,  =  0.55,  agreement  with  experiment  is  good 

everywhere  except  near  the  shock,  as  shown  in  Figure  4.8a.  Again,  the  angle  of 
attack  used  in  the  computation  (8.34°)  is  the  corrected  value  obtained  by  Harris  from 
the  measured  value  (9.86°)  using  a  linear  analysis  for  wind  tunnel- wall  effects.  For 
this  case,  the  flow  has  a  supersonic  bubble  located  well  forward  on  the  airfoil  upper 
siuface  and  is  slightly  separated  at  the  foot  of  the  shock.  The  computation  predicts 
the  shock  location  slightly  aft  of  experiment,  but  it  does  show  (from  Figure  4.8b)  the 
small  shock-induced  separation  evident  in  Harris's  measurements  between  x/c  =  0.1 
and  0.2. 


Figure  4-6.  Grid  Distribution  (200  x  63)  for  a 
NACA  0012  Airfoil  for  Viscous  Transonic  Computation 
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9. 


present  solution 


(a)  Pressure  Coefficient 


Figure  4-7.  Results  of  Viscous  Transonic  Flow  Over  a  NACA  0012  Airfoil 
at  a  =  1.49°,  M  =  0.7,  and  Re  =  9  x  10^ 


(b)  Density  Contours 


Figure  4-8.  Results  of  Viscous  Transonic  Flow  Over  a  NACA  0012  Airfoil 
at  a  =  8.34,  =  0.55,  and  Re  =  9  x  10^ 


Figure  4-9a  is  a  plot  of  pressure  coefficient  for  the  NACA  0012  airfoil  at  a  =  2.26®, 
=  0.749.  Again,  the  angle  2.26®  is  obtained  from  the  measured  angle  of  attack 
(2.86®)  using  a  linear  wind-tunnel-wall  correction  procedure.  For  this  flow  field,  a 
shock  wave  exists  on  the  airfoil  upper  surface  at  about  x/c  =  0.5,  which  is  strong 
enough  to  cause  significant  boundary-layer  separation.  For  this  case,  like  many 
others^^'^®^'^®^  the  present  theory  with  the  Baldwin-Lomax  model  does  not  capture 
the  shock  location  accurately,  and  predicts  its  location  too  far  aft  of  experiment. 

At  a  Mach  number  of  0.7,  several  angles  of  attack  were  computed  to  produce  the  lift 
coefficient  Cl  versus  angle  of  attack  a  plot  of  Figure  4-10.  The  experimental  data  of 
Harris  and  the  data  with  wind-tunnel-wall  correction  are  shown.  Agreement  with 
the  corrected  data  of  Harris  is  very  good. 

4.2.2  RAE  2822  Airfoil 

Two  studies  were  completed  for  the  RAE  2822  airfoil,  corresponding  to  Case  1  and 
Case  6  from  the  experimental  study  of  Cook,  et  A  144  x  69  O-grid  around  the 
airfoil  was  generated  by  an  algebraic  method,  and  is  shown  in  Figure  4-11.  The  outer 
boundary  is  located  at  a  distance  of  12  chord  away.  This  airfoil  is  a  supercritical 
airfoil  with  a  significant  amount  of  aft  camber.  Unlike  NACA  0012,  only  one 
experiment^ is  available  and  the  experimental  accuracy  is  unknown.  Case  1  of 
Reference  104  was  simulated  at  =  0.676,  Re  =  5.7  x  10^,  and  a  corrected  angle-of- 
attack  of  1.93°.  The  surface  pressure  coefficients  for  this  subcritical  case  in 
comparison  to  the  experiment  are  given  in  Figrire  4-12.  As  can  be  seen,  both  results 
are  in  good  agreement. 
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(a)  Pressure  Coefficient 


(b)  Density  Contours 


Figure  4-9.  Results  of  Viscous  Transonic  Flow  Over  a  NACA  0012  Airfoil 
at  a  =  2.26,  =  0.749,  and  Re  =  9  x  10^ 
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A  present  prediction 

o  Harris  (Linear  Correction  for  Interference) 
^  Harris  (uncorrected  data) 


angie-of-attack,  a 


Figure  4-10.  Lift  Coefficients  (Cl)  vs.  Angle-of-Attack  for  Transonic  Flow 
Over  a  NACA  0012  at  =  0.7,  and  Re  =  9  x  10^ 
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1 


Figure  4-11.  O-Grid  (144  x  64)  Around  RAE  2822 
Airfoil  For  Viscous  Transonic  Computations 
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Pressure  Coefficient 


Figure  4-12.  RAE  2822  Airfoil  Surface  Pressure  Distribution  at 
=  0.749,  Re  =  9  X  10^,  and  a  =  1.93°  (Case  1) 
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Case  6  was  computed  at  =  0.725,  Re  =  6.5  x  10^,  and  a  corrected  angle-of-attack  of 
2.79®.  Figure  4-13  shows  the  surface  pressure  coefficients  in  comparison  to  the 
experiment.  The  present  method  agrees  fairly  well  everywhere,  although  the  shock 
location  is  predicted  slightly  forward  of  the  experiment. 

4.3  NACA  0015  Airfoil  Undergoing  Constant  Rate  Pitching  Motions 

In  this  case,  we  considered  a  NACA  0015  airfoil  pitching  about  a  fixed  axis  at  a 
constant  rate  from  zero  incidence  to  a  maximum  angle-of-attack  of  approximately 
60®.  This  particular  airfoil  section  was  selected  because  many  experimental  studies 
are  available.^®^^®®  To  accommodate  the  time-dependent  coordinate  transformation 
(moving  grid)  in  the  present  flow  simulation,  a  "rigid"  grid  attached  to  the  airfoil  is 
used.  The  advantage  of  this  approach  is  that  once  an  initial  grid  is  generated,  the 
physical  coordinates  (x,y)  and  grid  speed  can  be  easily  computed  from  the 

prescribed  airfoil  motion. 

The  free  stream  Mach  number  is  0.02  and  the  chord  Reynolds  number  is  4.5  x  10^ 
which  are  the  same  as  experiment  of  Walker,  et  The  laminar  flow  is  assumed, 

as  in  the  work  of  Rumsey  and  Anderson,^®  and  Visbal  and  Shang.^® 

This  laminar  assumption  can  be  justified  from  the  following  two  reasons.  First,  a 
suitable  model  for  transition  and  turbulence  is  not  currently  available  for  the 
present  complex  unsteady  flow.  Second,  for  the  high  pitch-rate  regime,  the  energetic 
forcing  motion  is  expected  to  temporarily  dominate  over  some  transition  and 
turbulence  effects.  The  favorable  agreement  shown  below  between  the  predicted 
and  experimental  flow  features  confirms  the  above  hypothesis. 
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Pressure  Coefficient 


Figure  4-13.  RAE  2822  Airfoil  Surface  Pressure  Distribution  at 
=  0.725,  Re  =  6.5  x  10^,  and  a  =  2.79°  (Case  6) 
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A  grid  of  200  x  77  with  local  view  shown  in  Figure  4-14  was  used  in  the  present 
simulation.  This  grid  size  is  compatible  with  the  grid  in  the  work  of  Rumsey  and 
Anderson  (257  x  97  C-grid)^®  and  Visbal  and  Shang  (203  x  101  O-grid)^®.  The  grid 
extends  15  chords  away  from  the  airfoil.  The  temporal  accuracy  study  has  also  been 
carried  out,  and  it  is  found  that  for  AtU^/c  =  0.002  with  first-order  backward  Euler  in 
time,  the  solution  is  time-step  independent.  This  time  step  was  used  in  all  the 
computations  reported  below. 

Figure  4-15  displays  the  instantaneous  streamlines  at  several  angles  of  attack  and  the 
corresponding  experimental  visualization  data  for  a  non-dimensional  pitch  rate  of 
k  =  ac/U^  =  0.2. 

Both  experimental  and  computational  data  show  a  separation  bubble  (vortex)  on  the 
leading  edge  of  the  airfoil  as  well  as  a  separated  flow  on  the  trailing  edge  at  an  angle 
of  attack  of  27°  (Figure  4-15a).  The  leading  edge  vortex  is  also  known  as  a  dynamic 
stall  vortex. With  the  further  increase  in  the  angle  of  attack,  Figure  4-15b,  the 
leading  edge  vortex  grows  in  size  and  converts  along  the  upper  surface.  At  the  same 
time,  a  secondary  bubble  appears  on  the  leading  edge.  At  a  =  59°,  the  dynamic  stall 
vortex  has  grown  to  a  size  comparable  with  the  airfoil  chord  length,  and  the  flow 
falls  into  deep  stall  regime  (Figure  4-15c).  During  the  whole  process,  the  flow  is  fully 
attached  along  the  airfoil  lower  surface.  In  general,  the  present  prediction  compares 
favorably  with  the  experiment  of  Walker,  et 

With  an  increase  in  the  pitching  rate  to  k  =  0.4,  the  generation  of  the  leading  edge 
vortex  is  delayed  until  an  angle  of  attack  of  42°  (Figure  4-16b).  The  trailing  edge 
separation  is  visible  at  a  =  28°  (Figure  4-16a).  The  vortical  flow  is  less  vigorous  at  a 
=  60°  (Figure  4-16c)  in  comparison  to  the  case  of  k=  0.2  of  Figure  4-15c.  Again,  the 
calculations  predict  well  the  experimental  phenomenon. 
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Figure  4-14.  LvK'ai  View  a  200  x  /  /  O-Grid  Around 
NACA  0015  :or  Constant  Pitch-Rate  Simulation 


Figure  4-15.  Comparison  of  Computed  Flow  Field  With  Experimeiu  ;or  NACA  0015 
Airfoil  at  Constant-Rate  Pitch,  Re  =  45,000,  k  =  0  1,  lOi'  ~~  pnd 
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It  should  be  emphasized  that  the  pressure-based  method  can  deal  with  highly 
compressible  and  incompressible  flows  with  equal  efficiency,  while  the  density- 
based  method  has  difficulties  for  the  simulation  of  low  Mach  munber  flows,  such  as 
the  experiment  in  Reference  106,  where  M  =  0.02.  It  is  due  to  this  fact  that  Rumsey 
and  Anderson,^®  and  Visbal  and  Shang^®  used  M  =  0.2  in  their  simulations. 


I  wTTonTXiJW  I  FuTCl  !itc3  iTfTTSvv  irTiwlTm^VTiT*  lufiTO 


Unsteady  calculations  were  performed  for  a  NACA  0012  airfoil  pitching 
harmonically  about  the  quarter  chord  with  the  following  relation  of  angle-of-attack. 


a  =  +  da  sin  xot  (4.1) 

a^  =  4.86°,  da  =  2.44° 

The  non-dimensional  reduced  frequency,  k,  defined  as: 


/t  = 


or 


(4.2) 


is  0.081.  The  free-stream  Reynolds  number  is  4.8  x  10^,  and  Mach  number  is  0.6. 
Both  200  X  79  and  100x69  O-grids  similar  to  Figure  4-6  are  used  in  the  present 
simulations.  Third-order  TVD  scheme  and  second-order  Crank-Nicolson  in  time 
are  applied.  The  time  step  is  10*^  second,  resulting  in  approximately  2,000  time  steps 
per  pitching  cycle.  The  Baldwin-Lomax  turbulence  model  is  applied  to  account  for 
the  turbulence  effects  at  this  high  Reynolds  number. 


Figure  4-17a  shows  lift  coefficients  as  a  function  of  angle-of-attack  for  two  mesh  sizes 
of  200  X  79  and  100  x  69.  The  curves  are  followed  in  a  counter-clockwise  sense,  i.e., 
increase  in  alpha  is  represented  on  the  lower  portion  of  the  plots.  Both  grids  give 
good  agreements  with  experimental  data  of  London. Another  computation  with 
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□  experiment 
— 100  X  69  grid,  At  =  10"^ 
-  -  200  X  79  grid,  At  =  10“^ 


angle  of  attack  a 
(a)  spatial  accuracy 


□  experiment 
—  200x79  grid.  At  =  10’^ 
-  200  X  79  grid,  At  =  5.0  X  10' 


angle  of  attack  a 
(b)  temporal  accuracy 


Figure  4-17.  Lift  Coefficient  for  a  NACA  0012  Airfoil  Undergoing 
Oscillatory  Motion,  Re  =  4.8  x  10^,  =  0.6,  and  k  =  0.081 


reduced  time  step  size  by  a  half  (At  =  5.0  •  10'®)  with  200  x  79  grid  (see  Figure  4-17b) 
overlaps  the  one  with  At  =  10'^,  which  indicates  the  time  step  independence  of  the 
present  solution. 

The  pressure  coefficients  at  several  angles-of-attack  during  a  pitching  cycle  are 
compared  with  experimental  measurement,^^®  as  shown  in  Figure  4-18.  The  T 
symbol  denotes  that  the  angle-of-attack  is  increasing,  and  i  denotes  decreasing.  As 
can  be  seen,  the  agreement  is  remarkable  considering  the  complexity  of  the  flow. 

The  generation  and  disappearance  of  a  shock  wave  on  the  airfoil  upper  surface  can 
be  identified  from  the  pressure  coefficients.  Figvne  4-19  plots  the  density  fields  at 
several  angles-of-attack  within  a  pitching  cycle.  During  the  oscillatory  motion,  there 
is  a  shock  wave  on  the  upper  surface  of  the  airfoil,  and  the  flow  over  the  lower 
surface  is  predominately  subcritical.  Both  pressure  distributions  and  density 
contours  indicate  that  the  shock  position  oscillates  over  approximately  25%  of  the 
chord. 


4.5  Parametric  Effect  on  Dynamic  Stall  on  the  Airfoils 

As  it  is  generally  known,  dynamic  stall  depends  on  many  parameters.  This 
subsection  intends  to  investigate  the  influence  of  some  of  these  parameters  and  to 
explore  the  capability  of  the  present  code  in  predicting  dynamic  stall  phenomenon 
under  various  conditions. 
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□  expQriment 


4.5.1  Compressibility  Effect 


Here  the  experiments  of  Carr  et  and  Ahmed  and  Chandrasekhara^*^  were 
simulated.  The  conditions  are  the  same  as  in  the  experiments: 

NACA0012  Airfoil  Chord  Length  c=0.0762m 

Mach  Number,  =  0.3, 0.4  Reynolds  Number  Re  =  5.4  x  10^ 

Oscillating  Frequency,  f  =  21.64Hz  Reduced  Frequency,  k=0.05 
Angle  of  Attack  as:  a  =  +  Aa  sin  (ot,  Ojj,  =  10®,  Aa  =  10® 

A  grid  containing  200  x  79  cells,  as  shown  in  Figure  4-6  is  used.  First,  the  steady  state 
conditions  are  simulated.  The  density  contours  from  the  experiments  of  Carr  et  al 
using  real  time  interferometry  and  from  the  present  calculations  are  shown  in 
Figure  4-20  and  Figure  4-21.  For  Figure  4-20,  =  0.4  and  a  =  0.0®,  whereas  for 

Figure  4-21,  M„=0.3  and  a  =  10.78°.  The  experimental  fringes  seen  are  the  constant 
density  contours.  The  stagnation  point  in  both  experiment  and  calculation  is 
characterized  by  the  convergence  of  circular  fringes  (contours)  which  appear  at  the 
leading  edge.  The  density  contours  are  symmetric  on  both  the  upper  and  lower 
surfaces  in  Figure  4-20,  which  is  appropriate  for  this  0°  angle  of  attack.  The  abrupt 
turning  of  the  density  contours  on  the  upper  surface  indicates  the  presence  of 
boundary  layer.  At  a  =  10.78°  in  Figure  4-21,  the  stagnation  point  moves  to  the 
lower  surface  of  the  airfoil.  The  concentration  of  fringes  and  contours  near  the 
leading  edge  indicates  strong  acceleration  in  that  region.  It  is  clearly  seen  that  the 
present  solutions  correctly  simulate  the  experimental  observations. 
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(b)  Present  Prediction 


Figure  4-20.  Comparison  of  Density  Contours  between  the  Present  Prediction 
and  Experiments  of  Carr  et  al.  for  Steady  Flow  over 
a  NACA0012  at  M„  =  0.4,  a  =  0.0,  and  Re  =  5.4  x  10^ 
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(a)  Experiment  of  Carr  et  al.  (AIAA-90-0017, 1991) 


Figure  4-21.  Comparison  of  Density  Contours  between  the  Present 
Prediction  and  Experiments  of  Carr  et  al.  for  Steady  Flow 
over  a  NACA0012  at  =  0.3,  a  =  10.78°,  Re  =  5.4  x  10^ 
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Next  the  dynamic  stall  condition  is  simulated  with  k=0.05  and  =  0.3.  The 
oscillatory  motion  is  solved  by  marching  in  time  with  a  constant  time  step.  The 
time  step  size  is  such  that  it  corresponds  to  1,200  time  steps  to  complete  one 
oscillating  cycle.  The  third-order  TVD  scheme  in  space  and  second-order  accurate 
scheme  in  time  are  applied.  The  development  of  the  unsteady  flow  field  is  shown 
in  Figure  4-22.  During  the  upstroke  (t),  at  a  =  4.3°  (Figure  4-22a)  the  flow  is  fully 
attached.  As  the  airfoil  reaches  the  proximity  of  the  static  stall  angle  ia  =  12.0°), 
(Figure  4-22b),  the  attached  flow  persists.  The  boundary  layer  on  the  upper  surface, 
however,  has  grown  considerably,  as  seen  from  the  density  contours.  At  a  =  20° 
(Figure  4-22e),  a  shallow  reversed  flow  region,  which  initiated  at  a  =  12.0° T  (Figure  4- 
22b)  at  the  trailing  edge,  expands  upstream  rapidly.  This  flow  reversal  grows  in  size 
and  propagates  upstream.  It  is  this  reversal  flow  that  introduces  an  abrupt 
separation  of  the  boundary  layers  and  the  subsequent  development  of  a  leading  edge 
vortex.  With  further  decrease  in  the  angle  of  attack,  the  flow  reattachment  process 
starts  from  the  leading  edge  downward  (see  Figure  4-22f-4-22k).  At  a  =  6.4°,  the  flow 
on  the  upper  surface  attaches  completely. 

The  local  amplified  view  of  density  contours  and  experimental  fringes  during  the 
reattachment  process  of  the  above  oscillation  cycle  are  shown  in  Figure  4-23.  The 
leading  edge  bubble  due  to  shear  layer  attachment  and  its  growth  can  be  observed. 

4.5.2  Effect  of  Mean  Angle  of  Attack 

The  effect  of  the  mean  angle  of  attack  was  studied  by  using  the  same  conditions  as 
the  above,  except  that  A  a  =  6°,  and  is  selected  as  6°  and  12°.  For  =  6°,  it 
corresponds  to  light  stall,  and  for  a  =  12°,  it  is  a  deep  stall.  The  computed  lift 
coefficients  as  a  function  of  a  in  an  oscillation  cycle  for  =  6°  and  =12°  are 
shown  in  Figure  2-24.  The  experimental  data  of  McCrosby^  are  also  given.  The 
comparison  is  fairly  good. 
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Density 


Streamlines 
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Density 


Streamlines 


(f)  a  =  19.8®i 


(g)  a  =  18.4°i 


(h)  a  =  14.7®i 


Figure  4-22.  (Continued) 
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Streamlines 


Figure  4-22.  (Continued) 


Figure  4-23.  Density  Contours  and  Experimental  Fringes  Around 
a  NACA0012  During  the  Reattachment  Process  ot  an  Oscillating  Cycle. 
M^=0.3,  k=0.05,  Re=5.4xl0“,  aj^=UV'  and  Aa  =  10° 
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5. 


THREE-DIMENSIONAL  STATIC  AND  DYNAMIC  STALLS  ON  THE  WINGS 


In  this  section,  numerical  results  for  flow  fields  and  tip  vortex  formation  are 
presented  for  several  wing  sections.  These  results  are  compared  with  the  available 
experimental  data.  The  different  flow  conditions  and  the  winj  geometries 
considered  are: 

1.  Rectangular  wings  in  subsonic  flow; 

2.  Static  and  dynamic  stalls  on  a  rectangular  wing; 

3.  Static  and  dynamic  stalls  on  a  forward  swept  wing; 

4.  Static  and  dynamic  stalls  on  a  swept  back  wing;  and 

5.  Static  and  dynamic  stalls  on  a  delta  wing. 

5.1  Steady  Flow  over  Rectangular  Wings 

In  order  to  validate  the  present  3-D  code,  and  the  coding  of  the  Baldwin-Lomax 
turbulence  model,  a  simple  non-lifting  case  is  presented  first.  This  case  involves  a 
subcritical  flow  (M^  =  0.5)  about  a  large-aspect-ratio  wing  composed  of  NACA0012 
airfoil  sections.  Because  of  the  large-aspect-ratio  characteristic,  the  symmetry  plane 
solution  at  both  the  wing  root  and  wing  tip  is  essentially  two-dimensional  and 
should  compare  favorably  with  the  two-dimensional  counterpart.  Such  a 
comparison  of  pressure  coefficients  from  both  2D  and  3D,  and  experimental  data  of 
Thibert  et  is  shown  in  Figure  5-1.  In  the  2D  calculation,  a  grid  of  100  x  40  is 
used.  For  3D,  it  is  100x40x2,  and  symmetric  boundary  conditions  are  applied  at  two 
lateral  planes.  The  agreement  for  this  easy  case  is  excellent. 
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□  Experiment 


- 3D  Calculation 

- 2D  Caculation 


Figure  5-1.  Pressure  Coefficient  Comparisons.  NACA0012  Airfoil  Section, 
Large- Aspect-Ratio,  (TR=1.0,  =  0.5,  a  =  0®,  and  Re  =  3x10^) 
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The  second  geometry  consists  of  a  rectangular  wing  with  a  square  tip  and  a 
NACA0015  cross-section  without  twist  or  taper.  The  aspect  ratio  (AR)  is  2.5.  The 
incoming  free-stream  has  a  Mach  number  of  M^=  0.16,  Reynolds  number  of 
Re=2.0xl0^,  and  an  angle  of  attack  of  a=ll®. 

A  single-block  O-H  grid  topology  is  used  for  this  calculation.  There  are  84  points  in 
the  periodic  direction,  16  in  the  span  wise  direction,  and  48  in  the  normal  direction. 
The  wing  is  represented  as  a  solid  block  containing  4  radial  (normal)  cells  and  14 
spanwise  cells.  The  surface  grid  and  grid  root  wing  are  shown  in  Figure  5-2a.  Its 
planform  is  given  in  Figure  5-2b.  There  are  four  cells  extending  out  the  wing  tip  to 
capture  tip  vortex. 

The  wing  root  is  taken  as  a  symmetric  plane  to  remove  the  need  of  refining  the  grid 
in  order  to  resolve  the  boundary  layer.  As  a  result,  the  flow  features  near  the  wing 
root  resemble  those  of  two-dimensional  flows. 

Figure  5-3  shows  the  computed  surface  pressure  distributions  at  several  spanwise 
stations  compared  with  the  experimental  data  of  Spivey  and  Moorhouse^^^'^^^.  The 
inboard  spanwise  stations  show  good  agreement  with  the  experiment  as  evident  in 
Figure  5-3b-d.  However,  on  the  suction  side  at  the  wing  tip.  Figure  5-3a,  the  pressure 
distribution  is  not  as  well  predicted.  This  may  be  due  to  the  details  of  the  tip-cap,  as 
discussed  by  Srinivasan  et 

Figure  5-4  shows  the  surface  oil-flow  pattern  on  the  upper  (suction)  and  .ower 
(pressure)  sides  of  the  wing.  This  surface  oil  flow  picture  is  generated  by  releasing 
fictitious  fluid  particles  at  one  grid  point  above  the  surface  and  by  restricting  these 
particles  to  that  plane.  The  three-dimensional  effect  is  evident  near  the  wing  tip. 
Figure  5-5  shows  the  three-dimensional  perspective  view  of  particle  traces  around 
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Figure  5-2a.  Pictorial  View  of  a  82*45*16  3-D  Grid 
around  a  Rectangular  NACA0015  Wing 


symmetry  plane 


trailing  edge 


free>boundary 


(a)  Three-Dimensional  Perspective  View 

Figure  5-2.  Rectangular  NACA0015  Wing  and  Grid  Distribution. 

(a)  Three-Dimensional  Perspective  View;  (b)  Planform  and  Surface  Grid 
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(c)  z/B  =  0.36  (d)z/B  =  0.84 


Figure  5-3.  Surface  Pressure  Distributions  for  Several  Spanwise  Stations 
and  Comparison  with  Experimental  Data  for  NACA0015  Rectangular 
Wing.  M„=0.16,  a=ll°  and  Re=2xl06 
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(b)  Pressure  side 

Figure  5-4.  Surface  Oil-Flow  Pattern  (Simulated)  for  the  Rectangular 
NACA0015  Wing.  M^=0.16,  a=ll°,  and  Re=2xl0^, 

(a)  Suction  side  and  (b)  Pressure  side 


magnitude 


the  wing  tip.  It  is  seen  that  the  tip  vortex  is  lifting  off  from  the  upper  surface  of  the 
wing  at  about  the  trailing  edge  position.  Figures  5-4  and  5-5  reveal  that  the  particles 
released  in  the  lower  surface  of  the  wing  cross  over  the  tip  region  into  the  low- 
pressure  region  of  the  upper  surface  of  the  wing.  These  particles  mix  with  the 
particles  released  on  the  upper  surface  and  together  they  define  a  tip  vortex  that  is 
distinct  from  the  rest  of  the  sheet. 

Figiue  5-6  shows  the  pressure  contours  on  the  wing  surface.  With  the  exception 
being  near  the  tip,  the  pressure  is  almost  spanwise-independent.  It  is  evident  that 
the  tip  vortex  also  causes  a  pressure  gradient  at  the  wing  tip. 

5.2  Static  and  Dynamic  Stalls  on  a  Rectangular  Wing 

The  operational  wings  are  three-dimensional,  and  are  likely  to  encounter  unsteady 
flow  in  three-dimensional  conditions.  This  subsection  studies  the  three- 
dimensionality  of  unsteady  separated  flows  on  a  simple  unswept  symmetric  wing. 
The  experimental  visualization  of  Adler  and  Luttges'*®,  and  Ashworth  et  will  be 
use  as  a  comparison  basis,  since  only  limited  experimental  data  are  available. 

The  same  flow  conditions  are  used  as  in  the  experiment^:  Re=42,000,  M=0.02,  and 
angle  of  attack  varies  harmonically  with  time: 

a  =  +  da  sin  (cot) 

with  (5.1) 


k  =  ^  =  0.93,  a„  =  15°,  da  =  10° 
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(a)  Suction  side 


Figure  5-6.  Surface  Pressure  Contours  for  the  Rectangular  NACA0015  Wing. 

M„=0.16,  a=ll°,  and  Re=2xl0^ 

(a)  Suction  side  and  (b)  Pressure  side 


The  computational  domain  contains  58x48x18  grids  in  the  V/  ^  directions, 

respectively.  Again,  there  are  four  cells  beyond  the  wing  tip  to  capture  the  tip 
vortex.  All  the  computations  reported  below  are  carried  out  with  3rd-order  TVD  in 
space  and  Crank-Nicolson  schemes  in  time. 

To  gain  confidence  on  the  present  code,  a  comparison  is  made  first  between  the 
particle  traces  of  the  present  prediction  and  experimental  smoke  visualization  of 
Reference  48,  for  static  stall  and  d)mamic  stall  at  a  =  18®,  as  shown  in  Figure  5-7.  For 
the  static  stall,  the  flow  is  fully  separated  from  the  upper  surface  of  the  airfoil  at  this 
spanwise  location  (0.98c  from  the  wing  tip).  This  is  well  predicted  by  the  present 
theory.  The  present  calculation  also  predicts  several  vortices  on  the  upper  surface, 
which  is  not  seen  in  the  experiment.  It  should  be  noted  that  these  vortices  appear  in 
the  figure  only  if  particles  are  released  in  that  region.  Experimental  visualizations 
are  obtained  by  introducing  smoke  traces  upstream  to  the  airfoil.  The  convected 
smoke  particles  may  not  travel  through  the  upper  wing  surface  so  that  the  details  of 
the  above  recirculating  flows  may  not  be  detected. 

Under  unsteady  conditions  (dynamic  stall  in  Figure  5-7)  there  is  a  leading  edge 
bubble  (or  dynamic  stall  vortex)  on  the  upper  surface  in  the  experiment.  The 
agreement  of  the  present  theory  with  experimental  visualization  is  remarkable.  To 
gain  a  three-dimensional  understanding  of  the  flow  field,  oil-flow  patterns  on  the 
upper  wing  surface  and  particle  traces  at  different  spanwise  locations  for  both  static 
and  dynamic  stalls  are  shown  in  Figure  5-8  and  Figure  5-9,  respectively.  From  the 
oil-flow  pattern  of  both  the  figures,  a  transition  from  flows  dominated  by  the  wing 
tip  vortex  to  those  dominated  by  the  leading  edge  vortex  is  observed.  The  flow 
beyond  Ic  inboard  is  essentially  two-dimensional,  which  is  characterized  by  the 
leading  edge  vortex,  secondary  vortex  and  territorial  vortices.  Near  the  wing  tip, 
however,  the  secondary  and  territorial  vortices  are  suppressed  by  the  tip  vortex. 


90 


A  B 

Static  Stall  Dynamic  Stall 

a  =  18°  a  =  18°i 

ct^  =  15°,  Aa  =  10°,  k  =  0.93 


(a)  Experimental  Visualization 


(b)  Present  Prediction  (Particle  Traces) 


Figure  5-7.  Static  and  Dynamic  Stalls  on  a  Rectangular  Wing  of  a  NACA0015 

at  0.98c  Inboard.  Re=42.000  and  M=0.02. 

(a)  Experimental  Visualization  (b)  Particle  Traces  of  the  Present  Prediction 
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Figure  5-8.  Particle  Traces  on  the  Upper  Surface  of  NACA0015 
Rectangular  Wing  and  at  Several  Spanwise  Locations  for  Static  Stall. 
Re=42,000,  M=0.02  and  a =18.3° 
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Wing  root 


Dynamic  Stall  a  =  18.3® 


Figure  5-9.  Particle  Traces  on  the  Upper  Surface  of  a  NACA0015  Rectangular 
Wing  and  at  Several  Spanwise  Locations  During  Dynamic  Stall. 
Re=42,000,  M^=0.02  and  a=18.3‘’  Downstroke  (i). 
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such  that  only  the  leading  edge  separation  is  visible  from  the  particle  traces.  The 
size  of  the  secondary  vortex  grows  from  the  wing  tip  toward  the  wing  root. 

To  illustrate  the  size  and  convective  characteristics  of  the  leading  edge  vortex. 
Figure  5-10  shows  the  particle  traces  at  0.67c  inbound  location,  from  the  present 
computation  (k=0.93).  The  experimental  visualization®^  at  k=1.0  are  also  shown  for 
the  comparison.  At  the  maximvim  angle  of  attack,  a=25®,  the  shear  layer  vorticity  at 
the  leading  edge  is  coalescing  into  distinct  vortex  patterns.  With  decreasing  alpha, 
(Figure  5-lOb,  a=18®;  Figure  5-lOc,  a=5®),  the  leading  edge  vortex  convects  toward  the 
trailing  edge  by  the  free  stream  flow.  With  the  continuation  of  pitching  motion 
(Figure  5-10d,e)  the  leading  edge  vortex  is  shed  from  the  trailing  edge  to  the  wake 
and  the  flow  reattaches  to  the  surface  of  the  wing. 

The  flow  patterns  at  an  even  further  inboard  location  of  1.5c  are  given  in  Figure  5-11 
for  the  present  computation.  The  leading  edge  separated  flow  at  a=23.3®  (Figure  5- 
llj)  breaks  up  into  two  bubbles  (Figure  5-lla)  convecting  downstream  along  the 
upper  surface  of  the  wing.  At  a=21.0®  downstroke,  (Figure  5-llb),  the  shear  layer  is 
separated  from  the  leading  edge.  As  the  angle  of  attack  decreases,  the  shear  layer 
starts  moving  towards  the  airfoil  upper  surface  (Figures  5-llb  and  5-1  Ic).  The 
reattachment  occurs  when  the  static  stall  angle  is  reached  (Figure  5-1  Id),  and  a 
bubble  is  formed.  The  shear  layer  further  separates  downstream.  Meanwhile,  the 
leading  edge  attachment  progresses  and  the  bubble  grows  in  the  size  and  convects  to 
the  trailing  edge  during  the  upstroke  portion  of  the  cycle  (Figure5-llf-j). 
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The  difference  of  flow  patterns  between  the  upstroke  and  dow  stroke  for  the  bubble 
development  are  due  to  the  hysteresis  effects  that  are  always  present  in  these 
unsteady  flows.  The  effect  of  a  tip  vortex  can  be  appreciated  by  comparing  Figures  5- 
10  and  5-11  in  the  upstroke  of  the  cycle.  Near  the  wing  tip  (Figure  5-10)  the  flow  is 
attached  during  upstroke  instead  of  separated  for  the  location  far  inboard  (Figure  5- 
11). 

Figure  5-12  shows  the  flow  patterns  and  pressure  contours  on  the  upper  surface  over 
a  pitching  cycle.  The  movement  of  the  separated  flow  region  is  clearly  seen  from 
the  oil-flow  pattern,  while  separation  and  reattachment  are  always  accompanied  by 
high  pressure  gradients. 

Figure  5-13  shows  the  lift  coefficient  verses  the  angle  of  attack  in  one  oscillatory 
cycle.  The  2D  coimterpart  is  also  given.  It  is  interesting  to  note  that  the  lift  curve 
slope  for  the  3D  wing  is  reduced  from  the  2D  value.  This  characteristic  has  also  been 
found  in  the  experiments'^^'^^.  It  is  noted  that  dynamic  stall  is  delayed  and  the 
amplitude  of  the  lift  drop  is  reduced  in  the  3D  wing  in  comparison  to  that  in  2D. 
This  delay  in  dynamic  stall  is  attributed  to  lower  effective  angles  of  attack  caused  by 
proximity  to  the  wing  tip.  During  the  downstroke,  the  suction  pressure  peak  is 
usually  associated  with  the  stall  vortex.  This  vortex  interacts  with  the  tip  vortex  and 
is  suppressed  near  the  wing  tip.  As  a  result,  the  pressure  on  the  upper  wing  surface 
is  relatively  high  compared  to  the  2D  counterpart  (see  Figure  5-12).  This  causes  a 
lower  lift  during  the  downstroke.  The  oil  flow  visualization  in  Figure  5-12  indicates 
that  the  tip  vortex  rolls  up  over  the  wing  in  a  roughly  30®  triangle  beginning  at  the 
tip  leading  edge. 
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Figure  5-12.  (Continued) 


•e  5-1 


2D  calculation 
3D  calculation 


angle  of  attack;,  a 

Figure  5-13.  Lift  Coefficient  versus  Angle  of  Attack  for  a  Rectangular 
NACA0015  Wing  Under  Oscillating  Motion.  Re  =  42,000,  M  =  0.02  and  k  =  0.93 
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Forward  swept  wings  are  found  to  perform  well  at  high  Mach  numbers  and  have 
aerodynamic  advantages  at  very  low  airspeeds  (laminar  They  are  much 

better  suited  as  laminar  flow  wings,  due  to  smaller  effective  sweep  angles  compared 
to  a  corresponding  backward  swept  wing.  Thus,  the  laminar  botmdary  of  a 
comparable  forward  swept  wing  is  more  stable  against  attachment  line  transition 
and  crossflow  instability.  This  section  focuses  on  the  three-dimensional 
characteristics  of  unsteady  flows  produced  about  a  swept  forward  wing  imder  static 
conditions  and  dynamic  pitching  conditions. 

The  flow  and  wing  geometry  are  such  that  they  simulate  those  in  the  experiment  of 
Ashworth  et  The  wing  has  a  NACA0015  cross-section,  and  30®  forward  sweep. 
The  tape  ratio  is  1.0  and  the  aspect  ratio  is  2.0.  The  boundary  conditions  and 
computational  grid  on  the  wing  surface  are  shown  in  Figure  5-14.  The  Reynolds 
number  is  Re=40,000  with  M=0.02.  An  O-H  grid  topology  is  applied  to  generate  the 
grid  around  the  wing,  as  shown  in  Figure  5-2.  There  are  58  grids  in  the 
circumferential,  44  in  the  radial,  and  18  in  the  spanwise  directions.  Four  cells  are 
located  beyond  the  wing  tip  to  capture  the  tip  vortex  (see  Figure  5-14). 
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wing  tip 


Figure  5-14.  Top  View  of  the  Computational  Grid  and  Boundary 
Conditions  on  a  30'’  Forward  Swept  Wing 
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5.3.1  Computations  for  Static  Conditions 


Figures  5-15  to  5-17  show  the  comparisons  of  flow  visualizations^^  with  the  present 
3D  static  calculations.  The  Reynolds  number  is  40,000,  and  the  angles  of  attack  are 
3°,  15®  and  27°  for  Figures  5-15,  5-16,  and  5-17,  respectively.  At  an  angle  of  attack  of 
3°,  the  flow  is  fully  attached,  and  there  is  only  a  slight  difference  in  the  flow 
characteristics  at  0.58c  inboard  and  at  1.15c  inboard  (Figure  5-15).  At  an  angle  of 
attack  of  15°,  leading  edge  separation  appears  in  the  suction  side  of  the  wing.  The 
size  of  the  separated  flow  varies  with  spanwise  locations.  For  example,  at  0.58c 
inboard  (Figure  5- 16a),  the  flow  is  reattached,  while  no  reattachment  is  visible  at 
1.15c  inboard  (Figure  5-16b).  At  even  higher  angles  of  attack,  (a  =  27°)  the  flow  is 
totally  separated  from  the  upper  surface  (Figure  5-17).  It  can  be  seen  that  the  present 
computations  predict  well  the  experimental  visualizations. 

To  appreciate  the  three-dimensionality  of  the  present  simulation,  particle  traces  are 
shown  in  Figures  5-18  and  5-19  for  a=27°  with  top  and  side  views.  The  top  view 
(Figure  5-18)  clearly  displays  the  effect  of  tip  vortex  and  three-dimensionality  of 
separated  flow  on  the  suction  side  of  the  wing.  From  the  side  view  (Figure  5-19)  one 
can  see  the  growth  of  the  separated  flow  region  from  the  wing  tip  toward  the  wing 
root. 

Given  in  Figure  5-20  are  the  oil-flow  pattern  on  the  upper  surface  of  the  wing. 
There  is  a  focus  point  at  1.5c  inboard.  In  comparison  to  Figure  5-18  of  the  three- 
dimensional  particle  trace,  we  see  that  the  flow  particles  released  near  the  wing  tip 
are  sucked  to  the  wing's  upper  surface,  and  take  off  at  the  focal  point  to  form  a 
helical  flow  pattern. 
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I’rcssure  Cui\loui  s  l‘arlicle  Traces  Experiment  Visualizaitons 


gure  5-15.  ITesenl  C  alculations  and  Eh)w  Visualizations  for  a  3D  Forward  Swept  Wing 

at  <x  =  3^  Ke  =  40,000 


Forward  Swept  Wing 


Particle  Traces  Colored  by  Velocity  Magnitude 


Figure  5-18.  Particle  Traces  on  Static  3-D  30  deg.  Forward  Swept  NACA0015  Wing 

Re  =  40,000  M=0.02  Alfa  =  27  deg. 


Particle  Traces  Colored  by  Velocity  Magnitude 


Figure  5-19.  Separation  on  Static  3-D  30  deg.  Forward  Swept  NACA0015  Wing 

Re  =  40,000  M  =  0.02  Alfa  =  27  deg.  _ 


Reattachment  line  of  leading 
edge  vortex 


Figure  5-20.  Oil  Flow  Pattern  on  a  30°  Forward  Swept  Wing 
at  Re=4xl0‘‘,  a  =  27°  and  M  =  0.02 


5.3^  Computations  for  Dynamic  Conditions 


The  simulations  were  further  made  for  dynamic  wing  undergoing  oscillatory 
pitching.  The  angle  of  attack  varies  as  equation  (5.1).  The  rotation  axis  is  the  line 
connecting  the  root  and  tip  1/4  chord.  Figiure  5-21  gives  the  oil  flow  pattern  on  the 
wing  surface  and  particle  traces  at  several  spanwise  locations  during  an  oscillatory 
cycle.  Unlike  symmetric  (rectangular)  wings,  a  forward  swept  wing  experiences 
pronovmced  effect  of  tip  vortex.  This  is  evident  from  the  large  area  occupied  by  the 
tip  vortex  at  a  =  25®  (Figure  5-21a)  and  at  a  =  22.3®  (Figure  5-21b).  On  the  other  hand, 
the  leading  edge  vortex  domain  consists  of  a  triangular  wedge  with  base  far  inboard 
and  apex  near  the  leading  edge  of  the  wing  tip.  During  the  downstroke,  the  tip 
vortex  penetrates  inboard,  and  prevents  the  convection  of  the  leading  edge  vortex 
into  the  wake  (see  Figure  5-21d).  Like  a  rectangular  wing,  there  is  a  rapid 
disappearance  of  the  leading  edge  vortex  at  a  =  5®.  During  the  upstroke,  the  tip 
vortex  breaks  down,  (Figure  5-21  f).  The  result  is  massive  helical  flow  on  the  wing 
surface  (Figure  5-21g).  At  the  mean  time,  the  broken  tip  vortex  grows  with  the  angle 
of  attack  (Figure  5-21h). 

It  is  seen  that  the  forward  swept  wing  produces  flow  fields  that  differ  significantly 
from  those  of  a  straight  symmetric  wing.  The  flow  near  the  wing  tip  delays  the 
separation  far  below  the  static  stall  angle  during  the  downstroke. 
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The  previous  studies  on  two-dimensional  airfoils  show  that  the  upward  motion 
generates  a  leading  edge  vortex  that  passes  over  the  chord  to  the  trailing  edge. 
Depending  on  the  reduced  frequency  k  value,  a  trailing  edge  vortex  is  elicited  with 
the  opposite  sign.  This  trailing  edge  vortex  often  causes  a  rapid  separation  of  the 
flow  from  the  airfoil  surface.  In  the  present  calculation,  there  is  little  evidence  of  a 
trailing  edge  vortex.  It  appears  that  the  tip  vortex  suppresses  the  initiation  of  the 
trailing  edge  vortex  (see  Figure  5-21);  therefore,  the  dramatic  flow  separation  does 
not  coincide  with  the  passage  of  the  leading  edge  vortex  into  the  wake.  The  present 
work  is  consistent  with  the  work  of  Gad-el-Hak^^,  Carta^®,  and  Ashworth  et  al.^^  in 
that  vortices  form  over  the  upper  surfaces  at  high  angles  of  attack  and  these  vortices 
simply  increase  or  decrease  in  size  as  pitching  is  introduced. 

5.4  Dynamic  Stall  on  a  Swept  Back  Wing 

Computations  were  made  for  the  same  conditions  as  those  for  the  forward  swept 
wing.  A  30°  swept  back  wing  was  analyzed  to  investigate  the  effect  of  sweep.  Both 
static  and  dynamic  stalls  were  calculated,  but  only  dynamic  stall  is  presented  below. 

Given  in  Figure  5-22  are  similar  to  those  graphics  in  Figure  5-21,  showing  the  oil- 
flow  pattern  at  the  upper  wing  surface  and  particle  traces  at  several  spanwise 
locations.  In  contrast  to  the  forward  swept  wing,  the  leading  edge  vortex  now 
becomes  large  along  the  spanwise  direction  toward  the  wing  tip.  The  interaction 
between  the  leading  edge  vortex  and  the  tip  vortex  still  exists  but  is  only  limited  to 
the  proximity  of  the  wing  tip.  The  generation  and  propagation  property  of  the 
leading  edge  vortices  is  similar  to  that  in  the  forward  swept  wing.  The  unique 
aspect  of  a  swept  back  wing  during  dynamic  stall  is  that  the  swept  back  generates 
flow  from  the  wing  root  toward  the  tip.  This  flow  motion  resists  the  inboard  flow 
produced  near  the  wing  tip  by  the  tip  vortex.  Vortidty  is  then  accumulated  into  a 
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wing  upper  surface  0.6c  inboard  1.1c  inboard  1 .62c  inboard 
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Figure  5-22.  Continued 
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Figure  5-22.  Continued 


structured  leading  edge  vortex  on  the  swept  back  wing  sturface  at  spanwise  points 
proximal  to  the  wing  tip. 

To  compare  the  characteristics  of  the  interaction  of  the  leading  edge  vortex  and  the 
tip  vortex  for  the  three  different  wings  (30®  swept  back  wing,  rectangular  wing,  and 
30®  forward  swept  wing),  the  surface  oil-flow  patterns  at  a  =  22®,  and  a  =  8®  during 
the  downstroke  are  shown  in  Figure  5-23  and  Figure  5-24.  It  is  seen  that  the  forward 
swept  wings  and  swept  back  wings  produce  spanwise  fluid  motion  along  the  top 
surface.  The  minimum  pressure  line,  corresponding  to  the  leading  edge  vortex  core 
line,  runs  roughly  parallel  to  the  wing  leading  edge.  Since  the  airflow  is  not 
orthogonal  to  the  line,  the  span-directional  velocities  are  produced.  The  influence 
of  tip  vortex  is  the  most  pronounced  and  penetrates  far  inboard  for  the  forward 
swept  wing,  and  is  limited  only  to  the  wing  tip  for  the  swept  back  wing. 

Starting  from  a  wing  root,  the  leading  edge  vortex  size  increases  toward  the  wing  tip 
for  the  swept  back  wing,  keeps  constant  for  the  rectangular  wing,  and  decreases  for  a 
forward  swept  wing.  A  swept  back  wing,  for  which  the  wing  surface  area  covered  by 
a  vortical  structure  is  significantly  larger  than  the  forward  swept  wing,  has  the 
potential  of  providing  the  ideal  fluid  dynamic  characteristics  for  unsteady  lift 
enhancement.  This  is  proven  in  Figure  5-25,  which  shows  the  lift  coefficients  for 
the  above  three  wings. 
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(a)  30®  Swept  Back  Wing 
a  =  21®i 


(b)  Rectangular  Wing 
a  =  21.2®i 


(c)  30®  Forward  Swept  Wing 
a  =  23.3®>l 


Figure  5-23.  Comparison  of  Oil-Flow  Patterns  (Simulated)  on  the  Upper  Surface  of  a 
30®  Swept  Back  Wing,  a  Rectangular  Wing,  and  a  30®  Forward  Swept 
Wing  During  Downstroke  of  an  Oscillation  Cycle  around  a  =  22° 
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(a)  30®  Swept  Back  Wing 
a  =  8.87®i 


(b)  Rectangular  Wing 
a  =  6.7®i 


(c)  30®  Forward  Swept  Wing 
a  =8.47®i 


Figme  5-24.  Comparison  of  Oil-Flow  Patterns  (Simulated)  on  the  Upper  Surfaces 
of  a  30®  Swept  Back  Wing,  a  Rectangular  Wing  and  a  30°  Foward  Swept  Wing 
During  Downstroke  of  an  Oscillation  Cycle  Around  a  =  8° 
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rectangular  wing 


/ 


angle  of  attack,  a 


Figure  5-25.  Lift  Coefficients  versus  a  for  a  Rectangular  Wing, 
a  30®  Swept  Back  Wing,  and  a  30°  Forward  Swept  Wing 
Undergoing  Oscillatory  Pitching 
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5.5 


In  a  steady  flow,  the  lift  of  a  two-dimensional  airfoil  is  contributed  mainly  by  the 
leading  edge  suction  peak.  The  lift  increases  with  increasing  angle  of  attack  imtil  the 
stall  angle  is  reached.  The  separation  on  the  upper  surface  will  then  reduce  the 
leading  edge  suction  peak  causing  the  lift  to  drop.  The  lift  producing  mechanism  of 
a  delta  wing  is  somewhat  different.  There  are  two  smooth  such  peaks  inward  of  the 
leading  edges.  These  peaks  are  produced  by  a  pair  of  stationary  leading-edge  vortices 
formed  by  separated  flow  on  the  low-pressure  side  of  the  wing.  Therefore,  the  lift 
on  a  delta  wing  is  created  by  the  separated  vortical  structures  rather  than  by  the 
attached  flow  over  a  convex  surface. 

In  this  section,  the  unsteady  flow  field  around  a  delta  wing  is  studied  by  numerical 
simulation.  The  delta  wing  has  a  45  degree  sweep.  The  Reynolds  number  based  on 
the  root  chord  is  1.7  x  10^.  Figure  5-26a  is  a  sketch  of  the  45®  delta  wing,  which  has  a 
NACA0012  profile  at  each  spanwise  section.  The  wing  is  pitched  around  the  quarter- 
chord  position.  The  angle  of  attack  varies  with  time  in  the  same  way  as  shown  in 
Equation  (5.1).  The  reduced  frequency  is  0.24.  The  computational  grid  is  50  x  39  x  18 
(circumferential  x  radial  x  spanwise).  The  upper  surface  grid  is  shown  in  Figure  5- 
26b. 
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The  oil-flow  patterns  for  static  delta  wing  are  shown  in  Figure  5-27.  At  an  angle  of 
attack  of  5**,  the  flow  is  attacked,  and  the  leading  edge  vortex  roll-up  can  be  identified 
from  the  separation  line.  This  roll  up  vortex  is  a  spiral  vortex  sheet  shed  from  the 
leading  edge  as  shown  in  Figiire  5-28.  At  a  high  angle  of  attack,  the  separation  on 
the  delta  wing  starts  from  the  leading  edge  corner.  This  is  different  from  a 
rectangular  wing  which  initiates  a  separation  at  the  trailing  edge.  The  separation 
zone  propagates  to  the  wing  root  with  increasing  a  (see  Figure  5-27c-d).  There  is  a 
focal  point  on  the  wing  surface  at  a  =  15®.  This  in  turn  generates  a  helical  flow 
leaving  the  wing  surface.  The  three-dimensional  particle  traces  for  this  angle  (a  = 
15®)  are  shown  in  Figures  5-29  and  5-30. 

Figure  5-31  shows  the  surface  oil-flow  patterns  for  the  delta  wing  at  several  angles  of 
attack  during  an  oscillatory  cycle.  At  the  early  stage  of  the  upstroke.  Figure  5-31a-b, 
there  is  a  tip  roll-up  characterized  by  the  separation  along  the  leading  edge.  This  tip 
roll-up  forms  a  spiral  vortex  sheet  as  shown  in  Figure  5-28.  With  an  increase  in  the 
angle  of  attack.  Figure  5-31  c,  the  separation  of  the  leading  vortex  starts  from  the 
trailing  edge  corner,  and  propagates  upstream  and  inward  (see  Figure  5-31c-e). 
During  the  downstroke,  the  flow  starts  to  reattach.  For  example,  at  a  =  16.5°i,  the 
flow  near  the  whole  leading  edge  (from  the  apex  to  the  comer  of  the  trailing  edge) 
becomes  reattached.  It  is  interesting  to  note  that  the  leading  edge  vortex  does  not 
convect,  rather  it  experiences  a  grow-decay  cycle. 


(c)a  =  15‘’ 


((1)0  =  20** 


Figure  5-27.  Surface  Oil-Flow  Patterns  (Simulated)  for  a  Static  Delta  Wing 
at  a  =  5°,  10°,  15°  and  20°,  Re  =  17,000,  and  =  0.02 
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PARTICLE  TRACES  COLORED  BV  UELOCITV  riAGNITUDE 


Figure  5-29.  Particle  Traces  around  a  Static  Delta  Wing 
M  =  0.02  Alfa  =  15  deg.  Re  =  17,000 


PARTICLE  TRACES  COLORED  BV  UELOCITV  flAGNITUDE 


Figure  5-30.  Particle  Traces  around  a  Static  Delta  Wing 
M  =  0.02  Alfa  =  15  deg.  Re  =  17,000 


(c)  a  =  19.2°? 


(d)  a  =  25°T 


Figure  5-31.  Oil-Flow  Pattern  (Simulation)  on  a  45°  Swept  Delta  Wing  Undergoing 
Oscillatory  Pitching.  Re  =  1.7  x  10“*,  k  =  0.24,  M  =  0.02,  =  15°,  and  Aa  =  10° 
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6. 


CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  Conclusions 

The  design  of  future  generation  combat  aircraft  for  increasing  level  of  performance 
requires  good  control  capability  at  the  high  pitch  rates  at  angles  of  attack  beyond 
maximum  static  lift.  This,  in  turn,  require  detailed  knowledge  and  exploitation  of 
the  highly  imsteady  vortical  flow  field  in  the  vicinity  of  the  vehicle.  The  present  six- 
month  SBIR  Phase  I  study  developed  several  new  methodologies  for  pressure-based 
Navier-Stokes  equation  solvers.  With  the  developed  techniques,  both  steady  and 
unsteady  separating  flows  are  analyzed  for  2-D  airfoil  and  3-D  rectangular,  forward 
swept,  swept  back  and  delta  wings.  Based  on  the  present  investigation,  the 
following  concluding  remarks  can  be  made. 

1.  The  pressure-based  method  was  demonstrated  to  be  efficient  (both  in 
terms  of  storage  and  computation  time  requirements)  for  flows  ranging 
from  subsonic,  transonic  to  supersonic,  and  from  incompressible  to 
compressible  flows. 

2.  The  presently  developed  TVD  scheme  for  convective  term  discretization 
requires  no  artificial  dissipation  and  can  properly  resolve  the  concentrated 
vortices  with  minimum  numerical  diffusion.  When  applied  to  transonic 
flows,  the  TVD  scheme  can  capture  shock  with  a  single  transition  point. 
The  property  was  demonstrated  for  the  density-based  methods  but  it  has 
never  been  shown  before  for  the  pressure-based  methods. 

3.  For  inviscid  supersonic  flow  over  a  NACA0012,  the  pressure-based  TVD 
method  is  as  accurate  as  density-based  TVD,  and  performs  better  than  the 
density-based  method  with  artificial  dissipation. 
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4.  For  steady  viscous  transonic  flow,  the  present  analysis  well  captures 
transition  points  and  recirculating  flows.  It  is  very  competitive  relative  to 
previous  density-based  analyses. 

5.  Comparisons  with  experimental  visualization  and  measurement  of  the 
present  solutions  for  2D  airfoil  undergoing  both  constant  pitch  rate  and 
oscillating  motion  are  very  good.  All  the  essential  physics  are  well 
preserved. 

6.  Validation  study  for  3D  steady  flow  proves  the  accuracy  of  the  present 
code. 

7.  The  interaction  between  leading  edge  vortex  and  tip  vortex  for  3D  forward 
swept  wing,  rectangular  wing,  swept  back  wing,  and  delta  wing  have  been 
studied. 

For  a  rectangular  wing,  the  wing  tip  vortex  dominates  the  outboard 
stations  and  interacts  with  the  leading  edge  vortex  at  nearly  right  angle. 

For  a  forward  swept  wing,  due  to  the  induction  of  spanwise  flow  toward 
the  wing  root,  the  effect  of  wing  tip  vortex  penetrates  deep  into  the  wing 
surface,  and  suppresses  the  convection  of  leading  edge  vortex.  In  addition, 
the  size  of  the  leading  edge  vortex  grows  from  the  wing  root. 

For  a  swept  back  wing,  the  tip  vortex  effect  is  only  limited  to  the  proximity 
of  the  tip. 
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For  the  delta  wing,  the  separation  of  the  leading  edge  vortex  starts  from 
the  comer  and  shifts  inboard  with  pitch,  but  never  truly  convects 
downstream  as  in  two  dimensions. 


6.2  Recommendations  for  Future  Work 

The  current  (Phase  I)  study  has  successfully  verified  the  capability  of  pressure-based 
method  in  predicting  steady  and  unsteady  flows  on  3D  aerospace  vehicles.  Further 
studies  are  recommended  to  refine  the  computational  model  and  to  investigate 
several  dynamic  stcill  control  schemes. 

1.  Refinement  of  the  Present  Code.  In  Phase  I,  due  to  limited  availability  of 
experimental  data  and  limited  project  time,  only  static  lift,  pressure  coefficients  and 
dynamic  visualizations  have  been  compared  and  validated  against  experiments.  In 
Phase  n,  systematic  comparisons  of  predicted  drag,  lift  and  moment  coefficients  will 
be  made  with  benchmark  experimental  measurements  which  are  now  undertaken 
at  NASA  Ames  Research  Center  for  3D  rectangular  wing  with  NACA0015  cross- 
section.  Since  turbulence  modeling  is  an  essential  part  of  the  simulation,  the 
assessment  of  Low  Reynolds  number  turbulence  model  of  / .  Y.  Chien,  and  Standard 
k-e  model  will  also  be  carried  out. 

2.  Dynamic  Stall  on  Double-Delta  Wing.  For  a  double-delta  wing  under  static 
condition,  there  exists  many  interesting  phenomena  as  shown  in  Figure  6-1.  These 
include  the  interaction  of  strake  vortices  and  wing  vortices,  asymmetry  of  vortical 
patterns,  vortex  bursting,  and  vortex  sheet  tearing”^*^^^.  Further  study  into  dynamic 
condition  during  maneuvering  will  shed  light  on  the  fluid  physics  and  will  provide 
better  controlling  techniques. 
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3.  Study  of  Dynamic  Stall  Control  Schemes.  With  an  understanding  of  basic 
fluid  physics,  and  the  well  validated  code,  the  various  controlling  concepts  will  be 
investigated  by  computer  simulations.  Several  novel  schemes  proposed  by  Young 
are  shown  in  Figure  6-2^^®. 

In  Phase  n,  some  of  the  recently  proposed  osncepts  will  be  considered: 

a.  Vortex  flap  concept.  Figure  6-3a.^^' 

b.  Apex  fence  flaps.  Figure  6-3b^^^.  These  devices  are  deployed  at  an  angle  to 
slender  delta  wing.  They  alter  the  vortical  flow  field  and  produce  an 
intense  suction  at  the  apex  which  enhances  the  lift  and  gives  a  nose  up 
pitching  moment.  At  high  angles  of  attack,  they  reduce  apex  lift  and 
produce  a  desirable  nose-down  pitching  momentum. 

c  Forebody  strake.  Figure  6-3d^^^.  These  strakes  are  conformally  stored  in 
the  forebody,  and  when  deployed,  force  asymmetric  vortex  shedding  from 
the  forebody,  generating  a  controlled  yawing  moment. 

d.  Spanwise  blowing^^^'  Figure  6-3e.  With  realistic  blowing  rates,  the  jet 
momentum  can  stabilize  the  leading  edge  vortices  and  produce  significant 
lift  increments  at  high  angles  of  attack. 

4.  Computational  Flow  Visualization.  Dynamic  stall  numerical  simulation 
creates  large  data  sets  which  are  difficult  to  analyze  with  existing  graphic 
postprocessing  tools  such  as  PLOT3D,  FAST,  EXPLOT,  FIELD- VIEW  etc.  Two  types 
of  graphic  tools  need  to  be  developed  to  process  and  validate  the  computational 
results: 
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WINO  tCAOINO  EDGE 
PL  APS 


f—mih  vortex  flap 


r'. 


REDUCED  LIFT 
AT  high  alpha 


(a)  vortex  flap  concept 


(b)  apex  fence  concept 


(c)  deployable  strake  concept  (d)  spanwise  blowing  concept 

Figure  6-3.  Potential  Vortex  Control  Concepts  for  Phase  n  Study 
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a.  Dynamic  image  generation  with  animation  capabilities,  capable  of 
generating  shadowgraphs,  Schlierens,  and  interferometry  images,  smoke 
traces,  etc.  Automatic  detection  and  display  of  critical  point  lines  and 
surfaces  (separation  lines,  recirculation  bubbles,  etc)  are  also  essential. 

b.  Graphical  image  examiner  for  alignment  and  comparison  of 
computational  versus  experimental  amd  computational  versus 
computational  flow  images. 

Both  of  the  above  packages,  CFD-VIEW  and  CFD-Image  are  currently  being 
developed  at  CFDRC  on  SGI  graphic  stations  and  will  be  adapted  for  the  proposed 
dynamic  stall  flow  analysis  study. 

The  Phase  n  study  will  produce  a  validated  3D  CFD  code  which  will  be  of  significant 
value  to  the  U.S.  Air  Force,  Federal  Aviation  Agency  and  aircraft  manufacturers. 
This  will  also  provide  a  strong  foundation  for  futher  research  and  development  of 
various  dynamic  stall  control  schemes  for  advanced  combat  aircraft. 
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